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Multi-objective optimization of high speed 
vehicle-passenger catamaran by genetic 
algorithm 


Part II 
Computational simulations 


Zbigniew Sekulski, Ph. D. 
West Pomeranian University of Technology in Szczecin 


ABSTRACT 


Real ship structural design problems are usually characterized by presence of many 
conflicting objectives. Simultaneously, a complete definition of the optimum structural 
design requires a formulation of size-topology-shape-material optimization task unifying 
the optimization problems of the four areas and giving an effective solution of the problem. 
So far, a significant progress towards the solution of the problem has not been obtained. 
An objective of the present paper was to develop an evolutionary algorithm for multi- 
objective optimization of structural elements of large spatial sections of ships. Selected 
elements of the multi-criteria optimization theory have been presented in detail. Methods for solution of the 
multi-criteria optimization problems have been discussed with the focus on the evolutionary optimization 
algorithms. In the paper an evolutionary algorithm where selection takes place based on the aggregated 
objective function combined with domination attributes as well as distance to the asymptotic solution, is 
proposed and applied to solve the problem of optimizing structural elements with respect to their weight 
and surface area on a high speed vehicle-passenger catamaran structure, with several design variables, 
such as plate thickness, scantlings of longitudinal stiffeners and transverse frames, and spacing between 
longitudinal and transversal members, taken into account. Details of the computational models were at 
the level typical for conceptual design. Scantlings were analyzed by using selected rules of a classification 
society. The results of numerical experiments with the use of the developed algorithm, are presented. They 
show that the proposed genetic algorithm can be an efficient tool for multi-objective optimization of ship 
structures. 
The paper is published in three parts: Part I: Theoretical background on evolutionary multi-objective 
optimization, Part II: Computational investigations, and Part III: Analysis of the results. 


Keywords: ship structure; multi-objective optimization; evolutionary algorithm; 
genetic algorithm; Pareto domination, set of non-dominated solutions 


SEAGOING SHIP HULL STRUCTURE 
MODEL FOR MULTI-OBJECTIVE 
OPTIMIZATION 


General 


Effectiveness of the developed evolutionary algorithm for 
the multi-objective optimization of seagoing ship structures 
has been verified by solving the multi-objective optimization 
problem for the midship segment of the passenger-car 
catamaran ferry, based on the Austal Auto Express 82 design 
developed by Austal [HANSA (1997)], [Significant Ships 
(1997)], Figure 14. Models developed for the multi-objective 
optimization are: (1) ship structural model, (2) optimization 


/ Fig. 14. The Auto Express 82 high speed vehicle-passenger catamaran 
model and (3) genetic model. "Boomerang " 
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Fig. 15. Main particulars of the Auto Express 82 high -speed vehicle-passenger catamaran 


Structural model of the ship hull structure 


Main particulars of the Austal Auto Express 82 vessel are 
given in Fig. 15. The general arrangement of the ship and her 
corresponding cross- and longitudinal sections are shown in 
Fig. 16. For seagoing ships structural design in its initial stage 
concerns the cylindrical and prismatic zone amidships. For this 
reason the analysis of the midship block-section (17.5 x 23.0 
x 11.7 m) was assumed. Bulkheads form boundaries of the 
block in the longitudinal direction. In the block nine structural 
regions can be distinguished. All regions are longitudinally 
stiffened with stiffeners; their spacing being different in each 
structural region. The transverse web frame spacing is common 
for all the regions. Both types of spacing, i.e. of stiffeners and 
transverse frames, are considered the design variables. The 
transverse bulkheads were disregarded to minimize the number 
of design variables. 

The structural materials are aluminium alloys of the 
properties given in Tab. 1. The 5083-H111 aluminium alloys 
are used for plate elements while 6082-T6 aluminium alloys 
are used for bulb extrusions. The plate thicknesses and the bulb 
and T-bulb extruded stiffener sections are assumed according to 
the commercial standards and given in Tab. 2 ,3 and 4. The bulb 
extrusions are used as longitudinal stiffeners while the T-bulb 
extrusions are used as web frame profiles. Practically, the web 
frames are produced by welding the elements cut out of the 
metal sheets. Dimensions of the prefabricated T-bar elements 
are described by the four following design variables: web height 
and thickness as well as and flange breadth and thickness. In the 
case of extruded bulb a single variable is sufficient to identify 
the profile, its dimensions and geometric properties. It reduces 
the computational problem and accelerates analysis. 


Tab. 1. Properties of structural material — aluminium alloys 


125 (for 5083-H111 alloy) [N/mm?] 
250 (for 6082-T6 alloy) [N/mm?] 


70,000 N/mm? 
0.33 
26.1 kN/n? 


Property 


Yield stress Ro; 


Young modulus E 


Poisson ratio v 


Density p 


The strength criteria for calculation of plate thicknesses 
and section moduli of stiffeners and web frames are taken in 
accordance to the classification rules [UNITAS (1995)]. It was 
assumed that bottom, wet deck, outer side and superstructure 
I and II are subject to the pressure of water depending on the 
speed and the navigation region. The main deck was loaded 
by the weight of the trucks transmitted through the tires, the 
mezzanine deck the — weight of the cars, while the upper deck 
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— the weight of equipment and passengers. Values of pressure 
were calculated according to the procedures taken from the 
classification rules. 

In the study a minimum structural weight (volume of 
structure) and total outer area of structural elements intended for 
maintenance (cleaning, painting, etc.) were taken as the criteria 
and introduced to the definition of the objective functions 
and constraints defined on the basis of the classification 
rules. When structural weight and surface area are chosen 
as the objective functions , their values depend only on the 
geometrical properties of the structure (if structural material 
is fixed). The assumed optimization task is rather simple but 
the main objective of the study was to build the computational 
method, verify the computer code and prove its application to 
the multi-objective optimization of a ship hull. 


Tab. 2. Thickness of plates 


Thickness t [mm] 


OS} OO] A] HD] | |] UJN | = 


Tab. 3. Dimensions of bulb extrusions 


Cross-sectional area 
[cm?] 


Dimensions (h, b, s, s,)” 
[mm] 
80x19x5x7.5 
100x20.5x5x 7.5 
120x25x8x12 
140x27x 8x 12 
150x25x6x9 
160 x 29x 7x 10.5 
200 x 38x 10x 15 


D h-cross-section height; b - flange width; 
s - web thickness; s, - flange thickness. 


11.70 


9.25 


6.80 


23.00 


Fig. 16. Assumed model of craft structure — midship block-section, frame system and structural regions 
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Tab. 4. Dimensions of T-bulb extrusions 


Cross-sectional area 


Dimensions (h, b, s, s,)” 


[mm] 
200 x 100x 8x 15 
200 x 140x 8x 5 
200 x 60 x 10 x 12 
200 x 50x 8x 9.5 
210x 50x5x 16 
216 x 140x7.6x8 
220 x 80x5x8 
230 x 80x 10x8 
230 x 80x 5.8x8 
235x 170x 8x 10 
240 x 140x 6x 10 
260 x 90x 5x 9.5 
275 x 150x9 x 12 
280 x 100x 5x8 
280 x 100 x 8 x 10 
300 x 60x 15x 15 
310 x 100 x 7x 16 
310 x 123x5x8 
350 x 100 x 8x 10 
350 x 100x 5x8 
390 x 150x 6x8 
390 x 150 x 6x 12 
400 x 140x5x8 
410x 100x 6x 8 
420x 15x 5x10 
420x I$ x 8x 10 
450x 100x 9x 10 
450 x 150 x 10 x 12 


h — cross-section height, b - flange width, 
s - web thickness, s, - flange thickness. 


[cm?] 


2) 


Multi-objective optimization model of the ship 
hull structure 


Inthe most general formulation to solve ship structural multi- 
objective optimization problem means to find a combination of 
values of the vector of design variables x =[x, x, ... X, .. x,]" 
defining the structure which optimizes vector of the objective 
function f(x). The design variables should also meet complex 
set of constraints imposed on their values. The constraints 
formulate the set of feasible solutions. It is assumed that all 
functions of the multi-objective optimization problem are 
real and a number of constraints 1s finite. When considering 
computational costs an additional requirement may also be 
formulated that they should be as small as possible. 

For the multi-objective optimization problem a substitute 
scalar objective function may be formulated, on the basis of 
components of vector objective function, in the following 
form: 


F(x) = F(f (x), £(x)) = w fO) + w f(x) — min! (15) 


where: 
f(x) — astructural weight of midship block-section taken 


to optimization, 
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— anarea ofthe outer surface of structural members 
subjected to cleaning and painting operations 
(surface area for maintenance) in the section, 

w, and w, — weight coefficients used for partial optimization 

criteria. 


£(x) 


Taking into consideration the operational loads as well 
as the constraints imposed on the design variables especially 
those resulting from conditions of local and global strength 
formulated in the approved rules of a classification society 
the substitute scalar objective function can be expressed as 
augmented objective function of unconstrained minimization 
problem: 


f(x) = F(x) + 9 w,P, (x) = 
ae (16) 
=w f(x) + w,f(x)-* yw, Py (x) — min! 
k-l 
where all symbols are described before. 


The augmented objective function expression (Eq. 16) has 
been extended by components corresponding to dominance 
attributes and distance to the asymptotic solution. As 
a consequence the following form of combined objective 
function has been adopted: 


f(x) = w,u,(x) + wyu,(x) + 


W auc oa X) T LAE O) + (17) 
Ne 
+ W ossa 7 d.(x)] + Y wR (x) 
k=l 
where: 
w, and w, — weight coefficients used for partial optimization 
criteria, 
E — dominance rank coefficient, 
M — dominance count weight coefficient, 
"m — distance from asymptotic solution weight 
coefficient, 
u — utility function for structural weight: 
f, n -fa 
u(x) = fima HOY) — max! (18a) 
Ds 
u, — an utility function for area of the outer surface 
of structural members subjected to cleaning and 
painting operations: 
f, -ba 
u(x) = P ma - OD) — max! (18b) 
Dons 
where: 
f(x) — current value of the first optimization criterion, 
en — maximum value of the first criterion, 
£a) — current value of the second optimization 


criterion, 
— maximum value of the second criterion and all 
the remaining symbols are as outlined before. 


2,max 


As it has already been stated earlier, three aggregation-based 
multi-objective evolutionary strategies for taking account of 
the partial optimization criteria f (x) and f,(x) are used in the 
scalar objective function (Eq. 4) calculation, and therefore also 
in the fitness function value calculation (Eq. 16): 


— selection of variants by using the scalar objective function 
(Eq. 16) with the values of weight coefficients w, and w, 
set by the user (w_strategy = 2), 

— selection of variants by using the scalar objective function 
(Eq. 16) with the values of weight coefficients w, and w, 
randomly and independently generated in the range [0, 1] 
(w_ strategy = 4), 

— selection of variants by using the randomly selected single 
partial optimization criterion F(x) = w,f,(x) or F(x) = w,f,(x) 
(w_strategy = 3) which is implemented by the random 
selection of a single nonzero weight criterion. 
Additionally, it is also possible to have: 

— selection of variants without (w,,,, = 0) or with (W #0) 
(Eq. 17) by taking into account the dominance rank of 
feasible solutions, 


— selection of variants without (W un = 0) or with (Waun #9) 
(Eq. 17) by taking into account the dominance count of 
feasible solutions, 

— selection of variants without (W sancs — 0) OF with (Wisane £0) 
(Eq. 17) taking into account the distance of feasible 
solutions to the asymptotic solution. 


In the present formulation a set of 37 design variables 
is applied, cf. Tab. 5 and Fig. 17. Introduction of the design 
variable representing the number of transverse frames in 
a considered section: x , and numbers of longitudinal stiffeners 
in the regions: X,, X,, Xi» Xip X, X,,, X, X,,, X,,, enables to 
perform simultaneous optimization of both topology and 
scantlings within the topology-scantling optimization model. 


Tab. 5. Simplified specification of bit representation of design variables 


Substring length (no 


Value 


Lower limit Upper limit 


of bits) 


s X. 
imin imax 


1 serial No. of mezzanine deck plate 1 10 


- 


serial No. of mezzanine deck bulb 


serial No. of mezzanine deck T-bulb 


number of web frames 


number of mezzanine deck stiffeners 


serial No. of superstructure I plate 


serial No. of superstructure I bulb 


serial No. of superstructure I T-bulb 


BR pl l | l l | p 


number of superstructure I stiffeners 


serial No. of inner side plate 


serial No. of inner side bulb 


serial No. of inner side T-bulb 


number of inner side stiffeners 


serial No. of bottom plate 


serial No. of bottom bulb 


serial No. of bottom T-bulb 


number of bottom stiffeners 


serial No. of outer side plate 


serial No. of outer side bulb 


serial No. of outer side T-bulb 


number of outer side stiffeners 


serial No. of wet deck plate 


serial No. of wet deck bulb 


serial No. of wet deck T-bulb 


number of wet deck stiffeners 


serial No. of main deck plate 


serial No. of main deck bulb 


serial No. of main deck T-bulb 


number of main deck stiffeners 


serial No. of superstructure II plate 


serial No. of superstructure II bulb 


PT PS] PS] PST PS] PS] PST PS] PS] PS] PS] PST PS] PS] PST PS] PS] x xp Xx x xp x 


32 32 serial No. of superstructure II T-bulb 42 52 


number of superstructure II stiffeners 


serial No. of upper deck plate 


serial No. of upper deck bulb 


serial No. of upper deck T-bulb 


| | PT XL» 


number of upper deck stiffeners 


Multivariable string length 
(chromosome length) 


HHP W/ A SAIS HY] HY] AIS AIAI AIS AIAI AIS AIAI AISA S HPO! HW; A Ul AJAJ] Hwy 
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‘ Mezzianine deck 
' X,- plate thickness 
/ X, - longitudinal stiffeners bulb profile 
‘ x, - web frame T-bulb profile 
n x, - number of longitudinal stiffeners 


` Upper deck 

Xz,- plate thickness , 

Xə - longitudinal stiffeners bulb profile 

X, - web frame T-bulb profile 

Xy - number of longitudinal stiffeners 
^ 


Superstructure side Il 


Main deck 

X4, - plate thickness 

Xz - longitudinal stiffeners bulb profile 
X4, - Web frame T-bulb profile 

X» - number of longitudinal stiffeners 


Wet deck 
X» - plate thickness 

Xz - longitudinal stiffeners bulb profile 
X4, - web frame T-bulb profile 

X» - number of longitudinal stiffeners 


Inner side 
X, - plate thickness 


Xz - plate thickness 

X, - longitudinal stiffeners bulb profile 
Xs» - web frame T-bulb profile 

X3, - number of longitudinal stiffeners 


Superstructure side | 
x X,- plate thickness 
x, - longitudinal stiffeners bulb profile 
X, - web frame T-bulb profile 
X, - number of longitudinal stiffeners 


Transverse frame 
x, - number of web frames 


Outer side 

X4, - plate thickness 

X, - longitudinal stiffeners bulb profile 
X» - web frame T-bulb profile 

X4 - number of longitudinal stiffeners 


X, - longitudinal stiffeners bulb profile 


X,, - web frame T-bulb profile 


X,, - number of longitudinal stiffeners 


t 


Bottom 


X,, - plate thickness 

X,s - longitudinal stiffeners bulb profile 
X, - web frame T-bulb profile 

X, - number of longitudinal stiffeners 


Fig. 17. Assumed model of craft — specification of design variables 


Numbers of stiffeners and transverse web frames, 
varying throughout the processes of optimization, determine 
corresponding spacings. Scantlings and weights ofthe structural 
elements: plating, stiffeners and frames are directly depending 
on the stiffeners and frames spacings — topological properties 
of the structure. 

When optimizing the structural topology of the ship, 
a difficult dilemma is to be solved concerning a relation between 
the number of structural elements in longitudinal and transverse 
directions, and their dimensions, influencing the structural 
weight. And, should be also considered constraints related to 
the manufacturing process and functional requirements of the 
ship, e.g. transportation corridors, container supporting seats 
on the containerships (usually realized by longitudinal girders 
and floors in the double bottom) or positioning supports on 
the girders in the distance enabling entry of cars on ro-ro 
vessels. 

The behaviour constraints, ensuring that the designed 
structure is on the safe side, were formulated for each region 
according to the classification rules [UNITAS (1995)] 
constituting a part of the set of inequality constraints g,(x): 

* the required plate thicknesses t... , based on the permissible 
bending stress: 


t-te? (19) 
where: 
t = actual value of plate thickness in j-th region, 
* the required section moduli of stiffeners Z. 

Zj e LT = 0 (20) 
where: 
Z.,- actual value of the section modulus of stiffeners in j-th 


region, 
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e the required section moduli of web frames Z 


fj,rule* 


Z7 Leute Z 0 (21) 
where: 
Z,,— actual value of the section modulus of web frames in 
j-th region, 
* therequired shear areas of stiffeners Ade 
As E Pe 20 (22) 
where: 
Dm actual value of shear area of stiffeners in j-th 
region, 
* the required shear areas of web frames A, ṣe: 
A B ve aie z0 (23) 
where: 
Ag actual value of the shear area of web frames in j-th 
region. 


He side constraints h,(x), mathematically defined as 
equilibrium constraints, for design variables are given in Tab. 5. 
They correspond to the limitations of the range of the profile 
set. Some of them are formulated according to the author's 
experience in improving the calculation convergence. 

The additional geometrical constraints were introduced due 
to the *good practice" rules: 

* the assumed relation between the plate thickness and web 
frame thickness: 


t = Co 20 (24) 
where: 
t = actual value of the plate thickness in j-th region, 


[ue actual value of web frame thickness in j-th region, 


* the assumed relation between the plate thickness and 
stiffener web thickness: 


(25) 
where: 


t — 
j 
t .— 


s,wj 


actual value of the plate thickness in j-th region, 
actual value of stiffener web thickness in j-th region, 


e the assumed minimum distance between the edges of frame 
flanges: 


I(x,+1)- b; >0.3 m (26) 


where: 
b= 


- actual value of frame flange breadth in j-th region. 


The relationships supplement the set of inequality 

constraints g(x). 

] * B i$. dE 
After formulating constraints it is necessary to formulate 
mathematical form of the penalty function P,. The choice of 
mathematical form of penalty function is actually free, however 
certain basic requirements must be met: 

— promoting (preferring, awarding) the solutions which do 
not violate constraints, increasing in this manner the value 
of fitness function (selection probability) of the solutions, 

— penalizing the solutions which do not fulfil (violate) the 
constraints, decreasing , in this manner , the value of fitness 
function and consequently selection probability, 

— normalizing the value of fitness function to one. 


After conducting series of test computations, the identical 
mathematical exponential form of penalty function was 
assumed for all constraints, Eq. 9. Weight coefficients w, 
allow for implementing comparative, in relation to the others, 
meaning of a given constraint identified by coefficient k. 

Finally by taking into consideration all the specified 
assumptions, the evolutionary multi-objective optimization 
model can be written as follows: 

— find the vector of design variables x = [x, ... x, ... X,]^ Xp 
i= 1, ..., 37 as shown in Tab. 5, 

— optimize the combined objective function f(x) — min! given 
by Eq. 17, 

— subject to behaviour constraints given by Eq. 19+23, side 
constraints given in Ta. 5 and geometrical constraints given 
by Eq. 24226, build a set of equality h(x) and inequality 
g(x) constraints, 

— exponential forms of penalty functions for representing of 
constraints violation, Eq. 9. 


Genetic model of the ship hull structure 
General 


The ship structural multi-objective optimization problem 
described earlier contains a large number of discrete design 
variables and also a large number of constraints. In such 
a case the GA seems to be especially useful. Solving of the 
optimization problem by using GA requires to formulate an 
appropriate optimization model. The optimization model 
specified earlier was reformulated into an optimization model 
according to requirements of the GA, which was further used 
to develop suitable procedures and define search parameters 
to be used in the computer code. 

The genetic type model should cover: 

— specification of chromosome structure, 
— specification of fitness function: fitness — max!, 


— specification of genetic operators suitable for the defined 
chromosome structures and optimization task, 
— specification of the searching control parameters. 


Chromosome structure 


In this work, the space of possible solutions is the space 
of possible structural variants of the assumed model. The ship 
hull structural model is described (identified) as a collection of 
37 design variables, x,, described above. Each of them can be 
represented by a string of bits. For example the deck structural 
model is described (identified) as a collection of five design 
variables, x,, presented on Fig. 18. As we mentioned above 
the bit string is used as chromosome in GAs. Defining the 
chromosome structure consists in assuming: 

— a sequence in which variables in chromosome will be 
coded, 

— anumber of genes for recording every variable, 

— resolution capability of coding actual variable values, 

— the lowest and the highest values of the variables. 


solution trial — (...,  .., Xo * s Xo Mom) 


floating point vector 


(..., 10,..., 5 46 14 8,...) 5 (... 1101 ... 101 101101 0110 0100...) + 


binary coding 
chromosome — individual 


Fig. 18. Design variables and chromosome structure: main deck 


Coding sequence is free to choose and it does not affect 
algorithm computational characteristics, but this knowledge 
is essential for proper genotype decoding to decisive variables 
values - phenotypes. 

Number of genes in chromosome requires more attention. 
Generally, the higher gene number the greater demand for 
memory and increased computing time. Therefore number 
of genes cannot be too high. On the other hand, greater 
chromosome length allows for examination of objective space 
in greater detail. It cannot be too small otherwise mutation and 
crossover operators would not be able to function effectively 
on extremely short chromosomes and the resolution will be 
lower. The compromise is set by the user according to his 
own experience aiming at required resolution and acceptable 
time of computing. In this case memory size is of secondary 
significance. Bibliography review and the author’s own 
experience suggest that usually one variable is coded in 
a chromosome section containing 5 to 20 genes. 

Resolution capability determines sampling density of 
objective space. Greater resolution capability enables to 
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examine more points in objective space which increases 
probability of detecting interesting extremes. On the other 
hand, increased resolution capability results in increasing 
a number of computations and extends simulation time. The 
author’s opinion is that resolution capability should be set on 
the lowest acceptable level. In case the algorithm does not detect 
any interesting solutions the increasing of resolution capability 
should be considered. 

The user sets the lowest and highest values of variables on 
the basis of his own experience, according to the task being 
solved and providing appropriate convergence of algorithm 
while leaving discretion in exploration of objective space. 

A simplified specification for bit representation of all 
design variables is given in Tab. 5. The solution variant can 
be represented simply by a string of bits. 

The space of possible solutions is a space of structural 
variants of the assumed model. The hull structural model is 
identified by the vector x of 37 design variables, x. Each 
variable is represented by a string of bits used as chromosome 
substring in GA. The simple binary code was applied. Such 
coding implies that each variant of solution is represented 
by abit string named chromosome. Length of chromosome 
which represents of structural variant is equal to the sum of all 
substrings. Number of possible solutions is equal the product 
of values of all variables. In the present work the chromosome 
length is equal to 135 bits making the number of possible 
solutions equal approximately to 10?*. 


Fitness function 


A fitness function is used to determine how much the ship 
structure is suitable for a given condition in the optimum 
design with a GA. Because the combined objective function 
f(x) expressed by the relation Eq. 13 is: (1) well defined, (2) 
single-valued, (3) ascending, having real values and positive 
in the search space, it has been adopted directly to serve as the 
combined fitness function: 


f(x) = w,u,(x) + wu, (x) + 


W. aaO F Wooing) t 


count 


. Q7) 
T W distance 1 7 d.(x)] F EAS (x) 


k=1 


combined fitness = criteria + rank + count + 
+ distance + constraints 


All symbols used in the above given equation were already 
described. The expression (27) representing inclusion into 
selection process of test solutions, in addition to the degree 
with which they adjust to the established optimization criteria 
and formulated constraints, domination attributes as well 
as distance to the asymptotic solution, is the key point of 
the proposed Combined Fitness Multi-Objective Genetic 
Algorithm (CFMOGA). 


Genetic operators 


The basic genetic algorithm (Simple Genetic Algorithm 
- SGA) produces variants of the new population by using 
the three main operators which constitute the GA search 
mechanism: selection, mutation and crossover. In the present 
work the algorithm was extended by introduction of elitism 
and updating. 

The basic genetic algorithm (Simple Genetic Algorithm 
-SGA) produces variants of the new population by using 
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the three main operators which constitute the GA search 
mechanism: selection, mutation and crossover. In the present 
work the algorithm was extended by introduction of elitism 
and updating. Many authors described the selection operators 
responsible for chromosome selection due to their fitness 
function value [Goldberg and Deb (1991)], [De Jong (1995)], 
[Back (1996)], [Michalewicz (1996)]. After the analysis of 
the selection operators a roulette concept was applied for 
proportional selection. The roulette wheel selection is a process 
in which individual chromosomes (strings) are chosen according 
to their fitness function values; it means that strings with higher 
fitness value have higher probability of reproducing new strings 
in the next generation. In this selection strategy the greater 
fitness function value makes the individuals more important 
in a process of population growth and causes transmission of 
their genes to the next generations. 

The mutation operator which introduces random changes 
of the chromosome, was also described [Back (1996)], 
[Michalewicz (1996)]. Mutation is a random modification of 
the chromosome. It gives new information to the population 
and adds diversity to the mate pool (pool of parents selected for 
reproduction). Without the mutation it is hard to reach solution 
point which is located far from the current direction of search, 
while due to introduction of the random mutation operator the 
probability of reaching any point in the search space never 
equals zero. This operator also prevents against to the premature 
convergence of GA to one of the local optimum solutions, thus 
supporting exploration of the global search space. 

The crossover operator combines the features of two parent 
chromosomes to create new solutions. The crossover allows to 
explore a local area in the solution space. Analysis of the features 
of the described operators [Goldberg and Deb (1991)], [Back 
(1996)], [Michalewicz (1996)] led to developing a new, n-point, 
random crossover operator. The crossover parameters in this 
case are: the lowestn_x_site_minand the greatestn x site max 
number of the crossover points and the crossover probability 
p. The operator works automatically and independently for 
each pair being intersected (with probability p,), and it sets 
the number of crossover points n_x_site. The number of 
points is a random variable inside the set range [n x site min, 
n x site max]. The test calculations proved high effectiveness 
and quicker convergence ofthe algorithm in comparison witch 
algorithm realizing single-point crossover. It was also found 
that the number of crossover points n x site max greater than 
7 did not improve convergence of the algorithm. Therefore, the 
lowest and greatest values of the crossover points were set as 
follows: n x site min=1,n_x_ site max ^ 7. 

The effectiveness of the algorithm was improved with 
application of an additional updating operator as well as 
introduction of elitist strategy. 

Random character of selection, mutation and crossing 
operators can have the effect that these are not the best fitting 
variants of the parental population which will be selected for 
crossing. Even in the case they will be selected, the result 
will be that progeny may have a lower adaptation level. Thus 
the efficient genome can be lost. Elitist strategy mitigates the 
potential effects of loss of genetic material copying certain 
number of best adapted parental individuals to progeny 
generation. In the most cases the elitist strategy increases the 
rate of dominating population by well-adapted individuals, 
accelerating the convergence of the algorithm. The algorithm 
selects fixed number of parental individuals n, having the 
greatest values of the fitness function and the same number of 
descendant individuals having the least values of the fitness. 
Selected descendants are substituted by selected parents. 
In this way the operator increases exploitation of searching 


space. Update operator with fixed probability of updating p, 
introduces an individual, randomly selected from the parental 
population, to the progeny population, replacing a descendant 
less adapted individual. This operator enhances exploration 
of searching space at the expense of decreasing the search 
convergence. It also prevents the algorithm from converging 
toa local minimum. Both operators acts in opposite directions, 
and they should be well balanced; exploitation of attractive 
areas found in the searching spaces as well as exploration of 
searching space to find another attractive areas in the searching 
space depend on the user’s experience. 


Control parameters 


Single program run with the defined genetic model is 
characterized by values of eighteen control parameters. In this 
case the set of genetic model parameters set for each simulation 
run signed as symi includes 18 elements: 

symi = (n m no D, D, Pw Po c strategy, 
n x site min, n x site max, p,, elitist, 


W strategy, Wi Ww» Wank? W eoni W distance) 

where: 

ny — number of design variables (number of 
genes), 

P — chromosome length (number of bits), 

n, — number of generations, n, - size of 
population, 

n, — number of pretenders, 

Pin — mutation probability, 

P. — crossover probability, 

c_strategy — denotation of crossover strategy (0 for fixed, 


1 for random number of crossover points), 
n x site min — the lowest number of crossover points, 
n x site max — the greatest number of crossover points, 


D, — update probability, 

elitist — logical variable to switch on (elitism — yes) 
and off (elitism = no) the pretender selection 
strategy, 

w strategy — — denotation of strategy for aggregation of 


objective function, 

— weight coefficient of weight of structure, 

— weight coefficient of surface area of structural 
element intended for cleaning and painting, 

— weight coefficient of solution dominance 
rank, 

— weight coefficient of individual dominance 
count, 

— weight coefficient of distance of individual 
from asymptotical solution. 


WwW 
W 


1 
2 
rank 

count 


distance 


These 18 parameters control the successive simulation runs 
and identify them unambiguously for the adopted structure 
model. 

For selection of more control parameters it is not possible 
to formulate quantitative premises because of the lack of 
an appropriate mathematical model for analysis of GA 
convergence in relation to control parameters. The control 
parameters were set due to test calculations results to achieve 
a required algorithm convergence. 

All genetic parameters are specified by the user in advance 
of the calculations. This option is very important; the control 
of the parameter permits to perform search in the direction 
expected by the designer and in some cases it allows to find 
solution faster. The population size, number of variables and 
number of bits per variable, the total genome length, number 


of individuals in the population are limited by the available 
computer memory. 


Conclusion 


Finally, with taking into consideration all specified 
assumptions, the genetic model can be determined by the 
following: 

— chromosome structure specified in Tab. 5, 
— fitness function given by Eq. 28, 

— genetic operators, 

— control parameters. 


COMPUTATIONAL INVESTIGATIONS 
— THE SEARCH FOR SET OF NON- 
DOMINATED SOLUTIONS 


Computational investigations program 


In order to verify the suitability of the proposed method 
and the computer code developed for the seeking of Pareto- 
optimal solutions of the formulated multi-objective seagoing 
ship structure optimization problem, a number of calculation 
experiments have been performed, Tab. 7, by using the 
ship structure models earlier formulated and discussed in 
Section 6. 

From the multi-objective optimization point of view, the 
aim of the simulation was searching for non-dominated variants 
with respect to two optimization criteria with varying strategies 
for setting the values of weight coefficients for various criteria 
as well as dominance attributes: 

(1) Series 1.: the simulations marked with symbols sym1-1, 
sym1-2 and sym1-3 having the following values of control 
parameters: 

sym1-1: (37, 135, 10.000, 5.000, 10, 0.086, 0.800, 1, 1, 7, 
0.33, yes, 2, 0.5, 0.5, 0.0, 0.0, 0.0), 


sym1-2: (37, 135, 10.000, 5.000, 10, 0.086, 0.800, 1, 1, 7, 
0.33, yes, 3, random in [0,1], random in [0,1], 0.0, 0.0, 0.0), 


sym1-3: (37, 135, 10.000, 5.000, 10, 0.086, 0.800, 1, 1, 7, 
0.33, yes, 4, random 0 or 1, random 0 or 1, 0.0, 0.0, 0.0). 


In the simulation marked as sym1-1 the fixed values of 
weight coefficients are used for whole simulation: w, = 0.5 
and w, = 0.5; which refers to a classical method of weighted 
criteria. In the simulation marked as sym1-2 the values of 
weight coefficients w, and w, were generated by the computer 
code as random variables in the range [0, 1], which was done 
independently for each variant whenever the value of fitness 
function is calculated. In the simulation marked as sym1-3 the 
values of weight coefficients w, and w, were generated by the 
computer code as random variables equal to either 0 or 1, 
which was done independently for each variant whenever the 
value of fitness function is calculated; the value of 1 was used 
only for one, randomly selected criterion, with the remaining 
ones equal to 0. 

(2) Series 2.: the simulations marked with symbols sym2-1, 
sym2-2 and sym2-3 having the following values of control 
parameters: 

sym2-1: (37, 135, 10.000, 5.000, 10, 0.086, 0.800, 1, 1, 7, 
0.33, yes, 1, 0.0, 0.0, 3.0, 0.0, 0.0), 


sym2-2: (37, 135, 10.000, 5.000, 10, 0.086, 0.800, 1, 1, 7, 
0.33, yes, 1, 0.0, 0.0, 0.0, 3.0, 0.0), 


sym2-3: (37, 135, 10.000, 5.000, 10, 0.086, 0.800, 1, 1, 7, 
0.33, yes, 1, 0.0, 0.0, 0.0, 0.0, 3.0). 


Search for non-dominated variants while excluding the 
optimization criteria from the process of variant selection 
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Tab. 7. Control parameters of computational investigations 


Designation 


symi (n 


dv? 


w strategy, w,, W, W 


Specification 
he n, D D, Paw p, € strategy, n x site min,n x site max, p,, elitist, 


Msc 


rank? Mami 


(37, 135, 10,000, 5,000, 10, 0.086, 0.800, 1, 1, 7, 0.33, yes, 2, 0.5, 0.5, 0.0, 0.0, 0.0) 


(37, 135, 10,000, 5,000, 10, 0.086, 0.800, 1, 1, 7, 0.33, yes, 4, random in [0,1], random in [0,1], 0.0, 0.0, 0.0) 


(37, 135, 10,000, 5,000, 10, 0.086, 0.800, 1, 1, 7, 0.33, yes, 3, random 0 or 1, random 0 or 1, 0.0, 0.0, 0.0) 


(37, 135, 10,000, 5,000, 10, 0.086, 0.800, 1, 1, 7, 0.33, yes, 1, 0.0, 0.0, 3.0, 0.0, 0.0) 


(37, 135, 10,000, 5,000, 10, 0.086, 0.800, 1, 1, 7, 0.33, yes, 1, 0.0, 0.0, 0.0, 3.0, 0.0) 


w,—w,-—0.0(w strategy = 1) which was governed in particular 
simulations only by: (1) the value of the dominance rank of 
feasible solution, W,, — 3.0, Woount = 9-05 Waistance 7 9-9, in the 
simulation marked as sym2-1, (11) the value of feasible variant 
dominance count, w,,,,, = 3-0, Waan = 0.0, Waictance 7 9-0, in 
the simulation marked as sym2-2, (iii) the distance between 
the feasible variant and the asymptotic solution, w,....., = 3.0, 
Wank 7 9-95 Woount = 9-0, in the simulation marked as sym2-3. 
The purpose of the simulation series was to find out whether 
the developed tool is effective in case of evolution being 
governed only by (i) dominance rank, (ii) dominance count, or 
(iii) distance from a asymptotic solution. This refers to modern 
algorithms of evolutionary multi-objective optimization, where 
the evolution is governed only by dominance attributes. 
In all the simulations the functions of penalties imposed 
in the violations of constraints were active, w, #0, k = 1, 2, 
nyi, 
The computational investigations were carried out first of 
all for two-objective problems, as in this case it is possible to 
present obtained results graphically in a multitude of ways, 
which facilitate their interpretation and analysis. 


Results of computational investigations 


Results of computational investigations — Series1: 
the simulation marked with symbols sym1-1, 
sym1-2 and sym1-3 


The results of simulation sym1-1 are going to be discussed 
in the most detailed way by presenting the results characteristic 
1.00 
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(37, 135, 10,000, 5,000, 10, 0.086, 0.800, 1, 1, 7, 0.33, yes, 1, 0.0, 0.0, 0.0, 0.0, 3.0) 


for the developed method. For other simulations only the most 
important results will be presented. 

The Fig. 19 presents the evolution of macroscopic values 
characterizing the evolution of generated and evaluated 
population of ship structure solutions in the simulation sym1-1: 
(1) the greatest fitness function value f „p (2) the lowest 
distance of the feasible solution from the asymptotic solution. 
Multi-objective optimization of ship structure was performed 
with regard to the structure weight f, and the surface area for 
maintenance f, in case of fixed values of weight coefficients 
w, = w, = 0.5 used for optimization criteria. Dominance 
attributes: dominance rank, dominance count and distance 
from asymptotic solution, being excluded from selection. 
The figure shows a desired continuous rise of the greatest 
value of fitness function f „ indicating rising quality of the 
best generated test solutions. The highest values of fitness 
function saturate already in 1668 generation, which means 
that in the following generations no solutions were generated 
being better adapted in the sense of the fitness function used, 
and the computation resources were squandered replicating 
the non-dominated solutions set. The lowest distance between 
a non-dominated solution and the asymptotic solution changes 
during the evolution, but above the threshold of 87.38% of the 
highest value found during the simulation this takes place to 
a very small extent. For example, in the 857th generation the 
distance of the closest solution fg57(x) from the asymptotic 
solution is 1.114), see Fig. 20a, while in the successive 858th 
generation the distance of the closest solution fgss(x) to the 
asymptotic solution increased to 1.177, Fig. 20b. For the same 


sym1-1, 
w: = wa = 0.5 (w_strategy = 2), 
Wack = 0.0, Wow = 0.0, Was, = 0.0 


^ maximum fitness. 


—— minimal distance to asymptotic objective. 


r 


10000 


Generation 
Fig. 19. The results of multi-objective genetic optimization of ship structure with respect to structure weight f, and the area of element surface f, in case 
of fixed values of optimization criteria weight coefficients w, = w, = 0.5 with dominance attributes being excluded from selection (sym1-1); the curves 
present the evolution of a highest value of fitness function f „, the lowest value of non-dominated solution distance from a asymptotic one; the values are 
dimensionless and standardized in [0,1] range in relation to the highest values found during the simulation 


D All values of distance from the asymptotic solution are calculated in normalized space because it is impossible to calculate the Euclidean 
distance in physical objective space when different objectives are measured in different units. 
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Fig. 20. Change of a structure of non-dominated solutions set when moving from generation 857 (a) to generation 858 (b); during the genetic multi-objective 
optimization of ship structure with respect to structure weight f, and surface area f, in case of the fixed values of optimization criteria weight coefficients 
w, — w,— 0.5 with dominance attributes being excluded from selection (sym1-1); circles represent non-dominated solutions, red points represent non- 
dominated solutions closest from the asymptotic one; dimensionless values are normalized to the interval [0,1] in relation to the highest values in the ser? 7 
change of non-dominated solution set structure and change of non-dimensional (normalized) of distance of the nearest solution from asymptotic solution 
caused not by change of values of partial optimization criteria but only change of set of non-dominated solutions set structure can be observed 


solution the change of the distance to the asymptotic solution 
occurs only due to the change of the structure of the set of 
non-dominated solutions. On the other hand, for example, in 
the 2156th generation the distance of the closest solution f^ 56 
(x) from to the asymptotic solution is 1.187, Fig. 21a, while 
in the successive 2157th generation the distance of the closest 
solution f3157(x) from the asymptotic solution decreased to 
1.087, Fig. 21b. In this case the change of the set structure 
caused another solution, already present in the set, became 
the solution closest to the asymptotic solution. 

In Fig. 22 the evolution of the structure of non-dominated 
solutions set in simulation sym1-1 is shown by using the 
selected time-based cross-sections of this set, e.g. for 1, 
2000, 4000, 6000, 8000 and 10.000 generations as examples. 
A systematic growth of the size of non-dominated solutions set 
is apparent with 6, 6, 9, 12, 11 and 14 non-dominated solutions 
respectively in the consecutive time-based cross-sections as 
well as the desired evolution of this set in the direction of more 
advantageous values of partial optimization criteria. 


Fig. 23 presents a detailed structure of the non-dominated 
solution set of the last generation, presented in the physical space 
of objectives and normalized space of objectives. It can be seen 
that the set of non-dominated solutions including 14 variants 
of ship structure was found during the simulation. For each 
non-dominated variant the values of optimization criteria were 
specified as: f (x) — structural weight and f,(x) —cleaned/painted 
surface area. The designer may select for further development 
one of these variants or a group of them deemed by him to be 
the best. For the variant closest from the asymptotic solution 
which was found in 5116 generation, f5116(x), the distance 
equals 1.096 in the normalized objective space, the structural 
weight is f(x) = 1086.28 kN and the cleaned/painted surface 
area is f(x) = 7422.10 m’. This variant may be recommended 
if there is a need to select a single solution for the formulated 
ship structure multi-objective optimization problem: 


Fm = Es) = [ffs fzsuse9T = 


= [1086.28 7422.10] - [EN m?] 


» Normalization of the optimization objective values makes it possible to calculate the distance from asymptotic objective in the Euclidean 
sense in cases when the axes of the co-ordinate system represent the objectives denoted in various units. 
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Fig. 21. Change of a structure of non-dominated solutions set when moving from generation 2156 (a) to generation 2157 (b); during the genetic multi- 
objective optimization of ship structure with respect to structure weight f, and surface area f, in case of the fixed values of optimization criteria weight 
coefficients w, = w, = 0.5 with dominance attributes being excluded from selection (sym1-1); circles represent non-dominated solutions, red points represent 
non-dominated solutions closest to the asymptotic one; dimensionless values are normalized to the interval [0,1] in relation to the highest values in the set; 
change of non-dominated solution set structure and change of non-dimensional (normalized) of distance of the nearest solution from asymptotic solution 
caused not by change of values of partial optimization criteria but only change of set of non-dominated solutions set structure can be observed 
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Fig. 23. Detailed specification of non-dominated solutions set obtained during the genetic multi-objective optimization of ship structure with respect to 
structure weight f, and surface area f, in case of fixed values of optimization criteria weight coefficients w, — w, — 0.5 with dominance attributes being 
excluded from selection (sym1-1); black circles represent non-dominated solutions, red dots represent non- -dominated solutions closest from the asymptotic 
one; dimensionless values are normalized to the interval [0,1] in relation to the highest values in the set 


14 POLISH MARITIME RESEARCH, No 3/2011 


sym1-1: we w= 0.5 (w. strategy = 2), 


Surface area for maintenance, in m? 


Surface area for maintenance, in m? 


1000 1200 1400 1600 1800 2000 2200 2400 2600 1000 1200 1400 1600 1800 2000 2200 2400 2600 
Structural weight, in kN Structural weight, in kN 


sym-1: wi = w; = 0.5 (w. strategy =2), 
in 6000 generation, 
set of 12 non-dominated solutions. 


imitt * = w =05 aiii =2), 
in 4000 generation, 


(1281.80, 7203.27) 


XN ieMMMMD 


Surface area for maintenance, in m? 


Surface area for maintenance, in m? 


i 


1000 1200 1400 1600 1800 2000 2200 2400 2600 1000 1200 1400 1600 1800 2000 2200 2400 2600 


Structural weight, in kN Structural weight, in kN 


A. 


(1098 28.7422 11j 


Surface area for maintenance, in m? 
Surface area for maintenance, in m? 


1000 1200 1400 1600 1800 2000 2200 2400 2600 1000 1200 1400 1600 1800 2000 2200 2400 2600 
Structural weight, in kN Structural weight, in kN 


Fig. 22. History of the evolution of non-dominated solutions set during the genetic multi-objective optimization of ship structure with respect to structure 
weight f, and surface area f, in case of the fixed values of optimization criteria weight coefficients w, — w, — 0.5 with dominance attributes being excluded 
from selection (sym1-1); black circles represent non-dominated solutions, red dots represent non-dominated solutions closest from the asymptotic one 
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Fig. 24. History of creating a set of non-dominated solutions in the first six timeframes in which a change in structure of the set of non-dominated solutions 
has occurred during the genetic multi-objective optimization of ship structure with respect to structure weight f, and surface area f, in case of the fixed 
values of optimization criteria weight coefficients w, =w, = 0.5 with dominance attributes being excluded from selection (sym1-1); red circles represent 
non-dominated solutions, blue dots represent feasible solutions in present generations; in case of the first generation two existed already 
non-dominated solutions are coming from initial generation 
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Fig. 25. History of creating a set of non-dominated solutions in the last six timeframes in which a change in structure of the set of non-dominated solutions 
has occurred during the genetic multi-objective optimization of ship structure with respect to structure weight f, and surface area f, in case of the fixed values 
of optimization criteria weight coefficients w, = w, = 0.5 with dominance attributes being excluded from selection (sym1-1); red circles represent non- 
dominated solutions, blue dots represent feasible solutions in present generations 
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The author found it an interesting task to examine in which 
way the evolutionary optimization algorithm is builds a set 
of non-dominated solutions (approximation of Pareto set). Is 
the hypothesis correct, which states that during the evolution 
the set of feasible solutions as a whole shifts in direction of 
asymptotic solution, and is the set of non-dominated solutions 
approximating Pareto front is an edge of the set of feasible 
solutions? To verify the hypothesis, in Fig. 24 are presented 
feasible sets and approximation sets for time intersections 
corresponding to six first modifications of approximation set, 
i.e. 1, 2, 3, 5, 6 and 7 generations. In Fig. 25 are also presented 
approximation sets for time intersections corresponding six last 
modifications ofapproximation set, i.e. 5441, 5882, 7081, 8051, 
8712 and 9390 generations. Modification ofthe approximation 
set took place by including , to the successive approximation 
set, selected solutions, non-dominated in the ongoing 
feasible set and non-dominated by existing approximation set 
simultaneously. 

Figure 24 confirms, according to Fig. 18 that the algorithm 
builds the approximation set very intensively in early 
generations; in practice, in every following generation a new 
solution is being incorporated to the approximation set. The slow 
receding of the approximation set from the set of non-dominated 
solutions in the direction ofthe origin ofthe coordinate system 


Is visible. Set of feasible solutions maintains steady position 
in the objective space, preserving also similar range in this 
space. Fig. 25 confirms very slowly progressing building ofthe 
approximation set in final phases ofthe simulation; sequentially 
generated non-dominated variants are separated in spaces about 
1000 generations of evolution simulated by the algorithm. The 
approximation set is already separated from the set of feasible 
solutions and sequentially incorporated to its non-dominated 
solutions already overcoming significant distance from the 
feasible set and approximation set in the objective space. To 
overcome so the large distance in the objective space is possible 
in the way of very high attractive mutations in the chromosomes 
of individuals. Location of the set of feasible solutions and its 
range in the objective space do not undergo significant changes 
in course of evolution. Slightly larger concentration of set of 
feasible solutions in the later generations is visible presented in 
the Fig. 25, however confirmation if it is an actual or apparent 
effect requires further research. 

The presented results of the evolution of non-dominated 
solutions set indicate that non-dominated solutions set evolution 
can proceed also in such a way that the set of dominated feasible 
solutions preserves more less stable position in evaluation space, 
in turn discovered in sequence of non-dominated solutions are 
discovered sparse very good solutions found in course of 
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Fig. 26. Results of genetic multi-objective optimization of ship structure with respect to structure weight f, and surface area f, in case of random values 
of optimization criteria weight coefficients w, and w, in range of [0, 1] with dominance attributes being excluded from selection (sym1-2); a) the curves 


present the evolution of the highest value of, “fitness function jM 


max’ 


the lowest value of non-dominated solution distance from a asymptotic one, b) detailed 


specification of final non-dominated solutions set in objective space, c) detailed specification of final non-dominated solutions set in normalized objective 
space; black circles represent non-dominated solutions, red circles represent non-dominated solution closest from the asymptotic one; dimensionless values 
are normalized to the interval [0,1] in relation to the highest values in the set 
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Fig. 27. History of the evolution of the non-dominated solutions set during the genetic multi-objective optimization of ship structure with respect to structure 
weight f, and surface area f, in case of random values of optimization criteria weight coefficients w, and w, in range of [0, 1] with dominance attributes 
being excluded from selection (sym1-2); black circles represent non-dominated solutions, red dots represent non-dominated solutions closest from the 

asymptotic one 
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an algorithm simulated evolution; they create a set of non- 
dominated solutions receding in the act of the simulation from 
the set of feasible solutions to asymptotic ideal evaluations. This 
way of evolution the author suggests to call “drop by drop”, 
because it resembles rain drops falling on the ground; in this 
case the set of dominated feasible solutions represents a rain 
cloud while the set of non-dominated solutions, approximation 
set, represents drops that have fallen onto the ground. 

Figure 26a presents the evolution of macroscopic values 
characterizing the evolution of generated and evaluated ship 
structure variants population in simulation sym1-2 in case of 
random values of weight coefficients w, and w, in range of 
[0, 1] with dominance attributes being excluded from selection. 
The figure shows a desired continuous rise of the greatest value 
of fitness function f „ indicating rising quality of the best 
generated test variants. The highest values of fitness function 
saturate already in 1416 generation. The lowest distance 
between a non-dominated solution and the asymptotic solution 
changes during the evolution but above the threshold of 86,69% 
of the highest value found during the simulation is results in 
insignificant changes. 

In Fig. 27 the evolution of the structure of the non- 
dominated solutions set in simulation sym1-2 is shown by 
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using, as examples, the selected time-based cross-sections 
of this set, e.g. for 1, 2000, 4000, 6000, 8000 and 10,000 
generations. The size of non-dominated solutions set in the 
selected generations is apparent with 8, 9, 13, 12, 10 and 8 
non-dominated solutions and it does not change significantly 
in the course of simulation. 

Figures 26b and 26c presents the detailed structure of the 
non-dominated solutions set of the last generation produced in 
simulation sym1-2. It can be seen that a set of non-dominated 
solutions including 8 variants of ship structure has been 
found during the simulation. For the solution closest from 
the asymptotic solution, which was found in 6145 generation, 
£6145(X), the distance equals 1.088 in the normalized objective 
space, the structure weight is f (x)= 1113.66 kN and the surface 
area for maintenance is f(x) = 7361.45 m": 


f i27 féiasX) = [frg s) fies)" = 
= [1113.66 7361.45]" - [KN m?] 
Figure 28a presents the evolution of macroscopic values 
characterizing the evolution of generated and evaluated ship 
structure variants population in simulation sym1-3 in case of 


random values of weight coefficients w, and w, equal to 0 
or 1, used for optimization criteria with dominance attributes 


sym1-3. 
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Fig. 28. Results of genetic multi-objective optimization of ship structure with respect to structural weight f, and surface area f, in case of random values 
of weight coefficients w, and w, equal to 0 or 1 with dominance attributes being excluded from selection (sym1-3); a) the curves present the evolution 
of a highest value of fitness function f „, the lowest value of non-dominated solution distance from a asymptotic one, b) detailed specification of final 
non-dominated solutions set in objective space, c) detailed specification of final non-dominated solutions set in normalized objective space; black circles 
represent non-dominated solutions, red dots represent non-dominated solution closest to the asymptotic one; dimensionless values are normalized to the 
interval [0,1] in relation to the highest values in the set 
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Fig. 29. History of the evolution of non-dominated solutions set during the genetic multi-objective optimization of ship structure with respect to structure 
weight f, and surface area f, in case of random values of weight coefficients w, and w, equal to 0 or 1 with dominance attributes being excluded from 
selection (sym1-3); black circles represent non-dominated solutions, red dots represent non-dominated solutions closest from the asymptotic one 
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being excluded from selection. The figure shows a desired 
continuous rise of the greatest value of fitness function f... 
indicating rising quality of the best generated test solutions. The 
highest values of fitness function do not saturate before 8888 
generation. i.e. very close to the end of simulation. The lowest 
distance between a non-dominated solution and the asymptotic 
solution changes during the evolution, but above the threshold 
of 88.30% of the highest value found during the simulation it 
takes place only to a very small extent. 

In Fig. 29 the evolution of the structure of the non- 
dominated solutions set in simulation sym1-3 is shown by 
using, as examples, the selected time-based cross-sections 
of this set, e.g. for 1, 2000, 4000, 6000, 8000 and 10.000 
generations. Cardinality of non-dominated solutions set in 
the selected generations amounts to 6, 18, 12, 12, 11 and 15, 
respectively, and it does not change significantly in the course 
of simulation. 

Figures 28b and 28c present a detailed structure of a non- 
dominated solutions set of the last generation produced in 
simulation sym3. In can be seen that the set of non-dominated 
solutions including 15 ship structural variants has been 
found during the simulation. For the solution closest to the 
asymptotic solution which was found in 7611 generation, 
£761:(x), the distance of 1.123 in the normalized objective space, 
the structural weight is f(x) = 1153.68 kN and the surface area 
for maintenance is f(x) = 7381.57 m*: 


E nis B f, 00 B [frs f1)" = 
= [1153.68 7381.57]! - [kN m2] 


Results of computational investigations — Series2: 
the simulation marked with symbols sym2-1, 
sym2-2 and sym2-3 


Figure 30 presents the evolution of macroscopic values 
characterizing the evolution of generated and evaluated ship 
structure variants population in simulation sym2-1: (1) the 
greatest fitness value f... , and (2) the lowest distance between 
the feasible variant and the asymptotic solution. Multi-objective 
optimization of ship structure with respect to structural weight 
f, as well surface area for cleaning and painting, f,, in the 
case of the values of optimization criteria weight coefficients, 
w, =w, = 0.0, dominance count, W sun = 0.0 and the distance 
from asymptotic solution, w,.,,,., 7 0.0; evolution of 
generations is governed by dominance rank, w,,., = 3.0, and 
constraints components, w, # 0, k = 1, 2, ..., n. The figure 
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shows a desired continuous rise of the greatest value of fitness 
function f... indicating rising quality of the best generated test 
variants. There can also be seen a very early, reached already 
in 585 generation, saturation of the fitness function maximum 
value f... The smallest distance ofthe non-dominated solution 
from the asymptotic solution varies in course of evolution, but 
only in a small range over 79.42% of the greatest value found 
during the simulation. 

In Fig. 31 the evolution of the structure of the non- 
dominated solutions set in simulation sym2-1 is shown by using 
as examples the selected time-based cross-sections of this set, 
e.g. for 1, 2000, 4000, 6000, 8000 and 10.000 generations. 
Cardinality of non-dominated solutions set in the selected 
generations amounts to 4, 8, 6, 8, 8 and 10, respectively, and it 
does not change significantly in the course of simulation. 

Figure 32 presents the detailed structure of the non- 
dominated solutions set of the last generation, presented in 
a physical space of objectives and the normalized space of 
objectives. In can be seen that a set of non-dominated solutions 
including 10 variants of ship structure has been found during 
the simulation. For each non-dominated variant the values of 
optimization criteria have been specified as: f. (x) — structural 
weight and f (x) — cleaned/painted surface area. The designer 
may select, for further development, one of the variants or 
a group of them deemed by him to be the best. For the variant 
closest from the asymptotic solution, which was found already 
in 196 generation, the distance equals 1.064 in the normalized 
objective space, the structural weight is f(x) = 1105.95 kN and 
the cleaned/painted surface area is f(x) = 7345.11 m*: 


P us = fo) = [fr19s(X) f5:1990)]" = 
= [1105.95 7345.11] - [EN m] 


This variant may be recommended if there is a need to 
select a single solution for the formulated ship structure multi- 
objective optimization problem. 

Figure 33 presents the evolution of macroscopic values 
characterizing the evolution of generated and evaluated 
ship structure variants population in simulation sym2-2. 
Multi-objective optimization of ship structure with respect 
to the structural weight f, as well as surface for cleaning and 
painting f, in case of the values of optimization criteria weight 
coefficients, w, = w, — 0, dominance rank, Wax — 0.0 as well as 
the distance from asymptotic solution, W sistance — 0.0; evolution 
of generations is governed by dominance count, W sunt = 3.0, 
and constraints components, w, #0, Ic 2: a n. The figure 
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Fig. 30. The results of multi-objective genetic optimization of ship structure with respect to the structure weight f, and the area of element surface f, in case 
of the fixed values of weight coefficients w, = w, = 0.0, with dominance count and the distance from a asymptotic solution being excluded from selection: 
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count 
the highest value of fitness function f... and the lowest value of the non-dominated solution distance from a asymptotic one; the values are dimensionless and 
standardized in [0,1] range in relation to the highest values found during the simulation 
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Fig. 31. History of the evolution of the non-dominated solutions set during the genetic multi-objective optimization of ship structure with respect to 
the structure weight f, and surface area f, in case of the fixed values of weight coefficients w, = w, = 0.0, with dominance count and the distance from 
a asymptotic solution being excluded from selection: w „„ 7 0.0, w = 0.0, the evolution is governed by dominance rank, w „= 3.0, and constraints 
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components, (sym2-1); black circles represent non-dominated solutions, red dots represent non-dominated solutions closest to the asymptotic one 
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Fig. 32. Detailed specification of the non-dominated solutions set obtained during the genetic multi-objective optimization of ship structure with respect to 


the structure weight f, and surface area f, in case of exclusion from selection: optimization criteria, w, = w, = 0.0, dominance count, w 
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Fig. 33. The results of multi-objective genetic optimization of ship structure with respect to the structure weight f, and the area of element surface f, in case 


of the fixed values of weight coefficients w, = w, = 0.0, with dominance rank and the distance from a asymptotic solution being excluded from selection: w, 
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standardized in [0,1] range in relation to the highest values found during the simulation 


shows a desired continuous rise of the greatest value of fitness 
function f indicating rising quality of the best generated test 
variants. The figure also shows fairly late, reached in 7301 
generation, saturation of the maximum value of fitness function 
f ,, Minimum distance of the non-dominated solution from the 
asymptotic solution varies during the simulation but only in 
a small range over 83.18% of the greatest value found during 
the simulation. 

In Fig. 34 the evolution of the structure of the non- 
dominated solutions set in simulation sym2-2 is shown by 
using as examples the selected time-based cross-sections 
of this set, e.g. for 1, 2000, 4000, 6000, 8000 and 10.000 
generations. Cardinality of non-dominated solutions set in the 
selected generations amounts respectively to 4, 12, 12, 12, 13 
and 13, and it does not change significantly in the course of 
simulation. 

Figure 35 presents a detailed structure of the non-dominated 
solutions set of the last generation, presented in a physical 
space of objectives and the normalized space of objectives. In 
can be seen that a set of non-dominated solutions including 13 
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variants of ship structure has been found during the simulation. 
For each non-dominated variant the values of optimization 
criteria have been specified as: f (x) —structural weight and f,(x) 
— cleaned/painted surface area. The designer may select, for 
further development, one of these variants or a group of them 
deemed by him to be the best. For the variant closest from the 
asymptotic solution, which was found in 5533 generation, the 
distance equals 1.047 in the normalized objective space, the 
structural weight is f (x) = 1192.04 kN and the cleaned/painted 
surface area is f(x) = 7327.41 m*: 


227 £5530%) = [f 1553309) f2:553300]" = 
= [1192.04 7327.41]T- [kN nv] 


This variant may be recommended if there is a need to 
select a single solution for the formulated ship structure multi- 
objective optimization problem. 

Figure 36 presents the evolution of macroscopic values 
characterizing the evolution of generated and evaluated ship 
structure variants population in simulation sym2-3. Multi- 
objective optimization of ship structure with respect to the 
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Wrank = 0.0, Wocunr = 3.0, Wastanse = 0.0, 
in 4000th generation, 
set of 12 non-dominated solutions. 
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sym2-2, 
w = W; = 0,0 (w strategy = 1), 
W rank = 0.0, Wooune = 3.0, W astane = 0.0, 
in 2000th generation, 
set of 12 non-dominated solutions. 
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` sym2-2, 
W, = wg = 0,0 (w strategy = 1), 


W aes, = 0.0, Woaure = 3.0, Werstance = 0.0, 


in 10000th generation, 
set of 13 non-dominated solutions. 
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Fig. 34. History of the evolution of the non-dominated solutions set during the genetic multi-objective optimization of ship structure with respect to the 
structure weight f, and surface area f, in case of the fixed values of weight coefficients w, =w, = 0.0, with dominance rank and the distance from a asymptotic 


solution being excluded from selection: w 


rank 


=0.0,w 


distance 


= 0.0, the evolution is governed by dominance count, w 


count 


= 3.0, and constraints components, 


(sym2-2); black circles represent non-dominated solutions, red dots represent non-dominated solutions closest from the asymptotic one 
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Fig. 35. Detailed specification of the non-dominated solutions set obtained during the genetic multi-objective optimization of ship structure with respect to 


the structure weight f, and surface area f, in case of exclusion from selection: optimization criteria, w, = w, = 0.0, dominance rank, w 


from asymptotic solution, w 


distance 


= 0.0 and distance 


rank 


= 0.0; evolution is governed by dominance count, W „7 3.0, and penalty components, (sym2-2); black circles represent 


non-dominated solutions, red dots represent non-dominated solutions closest from the asymptotic one; dimensionless values are normalized to the interval 
[0,1] in relation to the highest values in the set 
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Fig. 36. The results of multi-objective genetic optimization of ship structure with respect to the structure weight f, and the area of element surface f, in case 


of the fixed values of weight coefficients w, = w, = 0.0, with dominance rank and the dominance count being excluded from selection: w 
= 3.0, and constraints components, (sym2-3); the curves present the evolution of 


the evolution is governed by the distance from asymptotic solution, w 


distance 


=0.0,w 


nk coun 


,7 0.0, 


the highest value of fitness function f... and the lowest value of non-dominated solution distance from an asymptotic one; the values are dimensionless and 
standardized in [0,1] range in relation to the highest values found during the simulation 


structural weight f, as well as surface area for cleaning and 
painting, f, in the weight coefficients values of optimization 
criteria w, = w, = 0, dominance rank, W, = 0.0 and the 
dominance count, Wun = 0.0; evolution of generations is 
governed by the distance from asymptotic solution, W istance 
= 3.0, and constraints components, w, # 0, k = 1, 2, ..., n. 
The figure shows a desired continuous rise of the greatest 
value of fitness function f „ indicating rising quality of the 
best generated test variants. The figure also shows fairly 
late, reached in 5714 generation, saturation of the maximum 
value of fitness function f „„ Minimum distance of the non- 
dominated solution from the asymptotic solution varies during 
the simulation but only in a small range over 83.25% of the 
greatest value found during the simulation. 

In Figure 37 the evolution of the structure of the non- 
dominated solutions set in simulation sym2-3 is shown by using, 
as examples, the selected time-based cross-sections of this set, 
e.g. for 1, 2000, 4000, 6000, 8000 and 10.000 generations. 
Cardinality of non-dominated solutions set in these selected 
generations amounts to 4, 10, 12, 11, 14 and 13, respectively, and 
it does not change significantly in the course of simulation. 


26 POLISH MARITIME RESEARCH, No 3/2011 


Figure 38 presents the detailed structure of a non-dominated 
solutions set of the last generation, presented in a physical 
space of objectives and the normalized space of objectives. In 
can be seen that a set of non-dominated solutions including 13 
variants of ship structure has been found during the simulation. 
For each non-dominated variant the values of optimization 
criteria have been specified as: f. (x) — structural weight and 
f(x) — cleaned/painted surface area. The designer may select, 
for further development, one of these variants or a group of 
them deemed by him to be the best. For the variant closest to 
the asymptotic solution, which was found in 5305 generation, 
the distance equals 1.085 in the normalized objective space, the 
structural weight is f (x) = 1060.03 kN and the cleaned/painted 
surface area is f(x) = 7485.93 m*: 


fona = £55030) = [fr:530503). £2:5395Q)T" = 


= [1060.03 7485.93]"- [kN m?] 


This variant may be recommended if there is a need to 
select a single solution for the formulated ship structure multi- 
objective optimization problem. 
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Fig. 37. History of the evolution of the non-dominated solutions set during the genetic multi-objective optimization of ship structure with respect to the 
structure weight f, and surface area f, in case of the fixed values of weight coefficients w, — w, — 0.0, with dominance rank and the dominance count being 
excluded from selection: w „= 0.0, w= 0.0, the evolution is governed by distance from asymptotic solution, w — 3.0, and constraints components, 


rank count. distance 


(sym2-3); black circles represent non-dominated solutions, red dots represent non-dominated solutions closest from the asymptotic one 
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Fig. 38. Detailed specification of the non-dominated solutions set obtained during the genetic multi-objective optimization of ship structure with respect 


to the structure weight f, and surface area f, in case of exclusion from selection: optimization criteria, w, = w, = 0.0, dominance rank, w 
= 0.0; evolution is governed by distance from the asymptotic solution, w 
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A method for evaluating theoretical 
and real operation of diesel engines 
in energy conversion formulation taking 
into account their operating indices 
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Gdansk University of Technology 


ABSTRACT 


The article proposes valuating the operation of an arbitrary diesel engine, based on the sample case of 
a ship s main engine in which energy conversion processes take place in a given time. The above operation 
is understood as the energy transfer to the screw propeller in the given time in which the energy conversion 
into work and/or heat and its further transmission take place. The here proposed method for evaluating 
the operation of the main engines installed in marine power plants consists in comparing the operation of 
these engines to a physical quantity the measuring unit of which is the joule-second (joule xsecond). A new 
term is introduced which bears the name of the theoretical engine operation and is the standard (ideal) 
operation which can be compared to the operation of real engines revealing different levels of wear. It was 
shown that the calculations of the theoretical operation defined in the above way cannot make direct use of 
commonly known theoretical Diesel and Sabathe cycles. Instead they should use the cycles modified by heat 
abstraction taking place in accordance with the isobaric, or isothermal process. Other new terms introduced 
in the article are: the degree of excellence of energy conversion to work, considered as the measure of 
excellence of engine operation, and the degree of engine operation dissipation, being the measure of its 
real operation. It is shown that if in time t of engine operation the case takes place that: L, = idem and L, 
= idem, then the engine operation dissipation is equal to its mechanical efficiency. 


Key words: operation, energy, diesel engine, ship's main engine 


INTRODUCTION 


The engine operation can be interpreted (in a descriptive 
formulation) as the result of the existence ofthe engine energy 
state which enables the engine to convert the energy E into 
heat and work in a given time t [2, 6, 9, 10, 13]. That means 
that when analysing the operation of these engines we have 
to take simultaneously into account the energy generated 
by the engines and the time of its conversion [6, 7]. In this 
interpretation the engine operation (in valuating formulation) 
can be compared to a physical quantity which is described 
by a numerical value and the measure unit bearing the name 
of the joule-second [joulexsecond]. The engine operation 
understood in the above way becomes worse when the engine 
wear increases. This means that the value of the operation of 
each engine will decrease with time, compared to the operation 
of the (ideal) reference engine. 

Obviously, the energy conversion into heat and work, which 
is observed in the working spaces of each internal combustion 
engine, may take place in different times. In practice, it is 
essential that the work done in a given time is possibly the 
largest, or the given work is done in the shortest possible 
time. In is also important in practice that the maximum heat is 
released during the combustion, while the lost heat is possibly 
the smallest. Ifthe most favourable energy conversion cannot be 


achieved, we conclude that the engine is in the partial capability 
state [3, 4, 10, 12]. Consequently, of high importance is the 
analysis of the engine operation and its comparison with the 
theoretical operation represented by the (ideal) reference engine 
operating in accordance with the theoretical cycle. Therefore 
there is a need for determining the theoretical operation of the 
ideal engine. An attempt to determine such theoretical operation 
will be presented using the marine engine as the sample case 
[6, 7], but it will also refer to other piston internal combustion 
engines. 


MAIN ENGINE OPERATION IN 
DETERMINISTIC FORMULATION 


In the deterministic formulation, the operation of the 
main engine can be analysed in two cases. In the first case 
we assume that the energy E transmitted by the engine to the 
screw propeller does not change in the time interval [0, t]. In 
that case we can assume that the engine operation (D) can be 
given by the formula: 


D=Et (1) 


In the second case, when the energy (E) transmitted by the 
engine to the screw propeller depends on time (t), the operation 
D(t) is to be calculated using the formula: 
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D(t)=D= [Ecod: Q) 


E —energy converted and transmitted during engine 
operation, 
[0, t] — time interval of engine operation. 


The energy is released in the form of work or heat. 
Therefore, according to formula (2) the engine operation can 
be interpreted as [4]: 

t t 
D, - [L@at ; D0 = [Q@adt (3) 
0 0 

where: 

D, — operation connected with the performed work (L), 

D, — operation connected with the heat transfer (Q), 

[t.t] -— time interval of engine operation, 

t— time, independent variable. 


The conditions of main engine operation are mainly 
determined by the external conditions in which the sea-going 
ship sails and tasks done by the ship users (the crew) [8]. 
These conditions and tasks are the reason why in different time 
intervals different amounts of energy are converted inside the 
working spaces of these engines. However, in each case the main 
engine operation is determined by their performance areas [6, 
9, 15]. These areas are limited by their velocity characteristics, 
such as the external and control characteristics. If these 
characteristics are mapped against the main engine propeller 
characteristics, which also belong to the velocity characteristics, 
then the ranges of operation of these engines are defined, which 
can be illustrated in the form of areas (Fig. 1). 

The real characteristics differ from their theoretical 
equivalents (Fig. 1) as the average torque (M_) of the piston 
internal combustion engine is the function of the average 
effective pressure (p,), i.e. 


M, — Cp, (4) 


and the effective pressure (p,) depends remarkably not only 
on the fuel charge, but also on the engine rotational speed (n) 
[2. 9, 13]: 


P. = RAG p Wa Po Ro To Lo M No No Nd © 


where: 


c, — coefficient determining constructional properties of 
the engine, 

AG, — fuel charge, 

W, — (net) calorific value of the fuel, 

p,  - pressure of the air supplied to engine working 
spaces, 

R  - universal (absolute) gas constant (R = R-M), 

H E s la á : 

T,  — temperature of the air supplied to the engine working 
spaces, 

L,  — theoretical amount of air needed for combustion of 1 
kg of fuel, 

X — excess air number, 

n,  -— efficiency of filling of the engine working space 
(cylinder), 

T — indicated efficiency, 

m, ~ mechanical efficiency of the engine. 


The main engine operation is time-limited in the overload 
area, especially in the torque overload area, the area 9'-3'-4-9-9* 
in the real characteristic (Fig. 1). The engine cannot also operate 
long in the area of pure speed overload, the area 6-7'-10'-9-6. 
Obviously, the engine operation is even more time-limited in 
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N= k;n" for WP,- idem 
WP,,47 idem 
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Fig. I. Speed characteristics of the real piston engine (main engine) 
installed on a sea-going ship [6]: 

* engine performance area (the area of permissible engine loads, 
the area of so-called stable engine operation): 1*-2*-9*-9-10*-10-8*-1*, 
* engine overload area: 2-3-3'-4-5-6-7'- 10'-9-2, 

* area illustrating the range of engine operation: 
1'-2'-9'-3'-4-5-6-7-10"-8'-T*, 

N, — effective power, N „— nominal (rated) effective power, n — rotational 
speed, h „~ maximal setting of the injection pump, h, — nominal (rated) 
setting of the injection pump, h, — continuous setting of the injection pump, 
h., — minimal setting of the injection pump, NR... — maximal setting 


min max 


of the controller, NR, — nominal (rated) setting of the controller, 
NR „ — minimal setting of the controller, WP... — most difficult sailing 


min max 


conditions (largest resistance of ship motion), WP, — nominal sailing 
conditions, WP... — easiest sailing conditions (smallest resistance of 


ship motion), N= kn" — real propeller characteristic illustrating the 
dependence of the power needed by the propeller on the rotational speed 
(n), k,— real power factor, m — coefficient depending on the ship s hull 
shape (for displacement hulls m = 3.0) 


the area of both torque and speed overload, the area 4-5-6-9-4 
in Fig. 1. 

The real performance area (working area) of the engine is 
limited by the curves: h s — idem, NR, = idem, hi = idem 
and NR,.. = idem. But the range of the main engine operation 
is smaller than the performance area, as the main engine can 
only operate in such a way that it secures delivering the energy 
needed by the screw propeller for executing a given task. 
The amount of the energy needed by the propeller in given 
operating conditions depends on the real propeller characteristic 
N= f(k, n), which is different for different sailing conditions 
WP. = idem. 

The WP conditions can change between the easiest (WP in = 
= idem) and the most difficult (WP... — idem). These conditions 
determine the range of main engine operation[5], the area: 
1'-27-9'-3'-4-5-6-7'-10'-8-l'in Fig. 1. 

The effective power N, generated by the main engine 
(Fig. 1) includes the information how fast the effective work 
L, will be done in the given time t [3, 9]. That means that the 
power characterises the engine operation in the aspect of the 
rate of energy conversion to work, taking into account various 
losses, in particular thermal loss. But the engine operation 
consisting in the energy conversion into work is not possible if 
earlier the energy has not been converted to heat in the engine 
working spaces [3, 4, 9, 14]. Therefore when analysing the 
main engine operation we should take into account the energy 
delivered to the engine in the fuel-air mixture, which is initially 
converted to heat (Q) and then to work (L), rather than pure 
engine power output. 


These considerations show that in case of the main engines 
installed on ships, as well as all other piston internal combustion 
engines, converting the chemical energy contained in the fuel- 
air mixture generated in the combustion chambers into heat, 
and then into the mechanical energy of the crankshaft-piston 
assembly makes it possible to generate the effective power N.. 
This power is to be generated in time t needed for executing 
a given task by the ship. That means that in order to execute 
this task in the given time, the effective work L -t is to be done 
by an arbitrary main engine. The realisation of this work is the 
effect of the torque (M.) generated on the crankshaft at a given 
rotational speed (n) of each piston internal combustion engine, 
including the main engine [2, 10, 15]. Therefore according to 
formula (2), the engine operation interpreted as the energy 
conversion into effective work L, in time t can be given as: 


t t 
D, = [L.@at =27 [nCOM, Cd (6) 
0 0 

The main engines used on container ships, bulk carriers, and 
tankers are two-stroke, low-speed engines, most frequently with 
uniflow scavenging, or lateral scavenging with inlet ports higher 
than outlet ports [1, 10, 14]. A sample indicator diagram recorded 
on such an engine is shown in Fig. 2a for the case of uniflow 
scavenging, and in Fig. 2b for the engine with lateral scavenging 
and inlet ports (h,) higher than outlet ports (h.), h; > h, . 

The presented real cycles show that the compression 
of the fresh charge (the mixture of the air delivered to the 
cylinder and the remnants of the exhaust gas which remained 
in the cylinder after the previous cycle) takes place after 
the cylinder scavenging process, consisting in the inflow of 
the fresh air with the simultaneous outflow of the exhaust 
gas, has ended. Until then, starting from the time instant 
of removing the exhaust gas from the cylinder, the heat is 
abstracted with the exhaust gas to the environment. Therefore 
theoretical cycles of diesel engines should take into account 
not only isochoric heat abstraction (like in the Diesel 
and Sabathe cycles), but also, depending on engine speed [7]: 
* heat abstraction at constant pressure (isobaric), in case of 

low-speed engines, 

* heat abstraction at constant temperature (isothermal), in 
case of medium- and high-speed engines. 


The above approach is possible if the assumptions are 
adopted that: 

* in low-speed engines the flow resistances change 
insignificantly and we can assume that the heat abstraction 
process takes place at the same pressure. 

* however, during the realisation of this process remarkable 
temperature change can be observed due to a relatively 
long time of heat transfer to the environment (compared 
to high-speed engines), before the abovementioned fresh 
charge compression is initiated in the cylinder. 


On the other hand, in medium-speed and, especially, high- 
speed engines such analyses are possible if we assume that: 

* the flow resistances change remarkably (compared to 
low-speed engines) and we cannot assume that the heat 
abstraction process takes place at constant pressure. 

e the temperature change observed during the realisation of this 
process is neglectable due to a relatively short time of heat 
transfer to the environment (compared to low-speed engines) 
before the fresh charge compression starts in the cylinder. 
Therefore we can assume in that case that after the isochoric 
heat abstraction is ended, the isothermal or quasi-isothermal 
heat abstraction process takes place in these engines before 
the fresh charge compression starts in the cylinder. 


DMP 


DMP 
Fig. 2. Working cycles of two-stroke engines with a) uniflow scavenging, 
b) lateral scavenging for h,>h,, 1 — closing the air inlet ports, 

5'— closing the exhaust gas exit valve (Fig. a), closing the exhaust gas 
outlet ports (Fig. b), 5 — opening the exhaust gas exit valve (Fig. a), opening 
the exhaust gas outlet ports (Fig. b), I — opening the air admission non- 
return valves, P- ambient pressure, V,- combustion chamber volume, 

V. — displacement volume, V — total cylinder volume, GMP — upper dead 
i centre of the piston, DMP — lower dead centre of the piston, 

p — pressure, V — volume 


To evaluate the operation of a real engine in given 
conditions, as described in Ref. [2], its operation should be 
compared with the theoretical operation. 


THEORETICAL OPERATION OF 
INTERNAL COMBUSTION ENGINES AND 
ITS PRACTICAL MEANING 


In case of diesel engines, of high importance is not only the 
possible and needed operation, which depends on the effective 
work L, and time (6), but also the theoretical operation, which 
makes it possible to assess how far the possible engine operation 
differs from the standard operation. To do this, we need to know 
the theoretical work (L,) and the time t of the operation. The 
work L, can be easily calculated based on the average theoretical 
pressure (p). The theoretical operation, in other words: the 
theoretically possible operation (D, ) can be determined using 
formula (2) and taking into account relation (6) as: 
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t 


D, = [L,@de (7) 


0 

The theoretical cycles of diesel engines include 
thermodynamic cycles of this type of engines, such as the 
Sabathe cycle and the Diesel cycle. They are used for comparing 
the ideal and real course of processes of thermal energy 
conversion into mechanical energy in arbitrary diesel engines 
(and for this reason are frequently referred to as the reference 
cycles). The theoretical Sabathe cycle is a relatively good 
description of the real cycles recorded in high-speed diesel 
engines, while the theoretical Diesel cycle well models the real 
cycles taking place in low-speed engines of this type [8]. 

The calculations of the theoretical work (L,) need the 
information of the work (L,,) done during one theoretical cycle 
of engine operation. Unfortunately, the theoretical Sabathe and 
Diesel cycles cannot be directly used in this case. Instead, we 
should use their modified versions which take into account the 
specifics of operation of particular internal combustion engine 
types. In case of two-stroke engines with uniflow scavenging, 
this specific results from the delay of the outlet vale closing, 
which makes the fresh charge compression start not after 
reaching DMP by the piston, but after the exhaust gas outlet 
valve has been closed (Fig. 2 a, point 5’). At the same time 
for engines with lateral scavenging and port heights h, > h,, 
this specific results from the delay in closing the air inlet ports 
by the piston with respect to the DMP (Fig. 2b, point 1). The 
theoretical cycles, (Fig. 3a and Fig. 3b), in which compression 
of the working medium also does not begin at the DMP, but 
later at point 1, reflect the specifics of the both types of engines. 
A separate problem is taking into account the fact that the 
exhaust gas decompression ends at point 5’ of the real diagram, 
and not at DMP (Fig. 2). 

The presented discussion shows that for low-speed engines 
the theoretical cycle in the form of the modified Diesel cycle 
additionally taking into account the isobaric heat abstraction 
can be used (Fig. 3a), while the cycle used for medium-speed 
engines is the modified theoretical Diesel cycle taking into 
account additional isothermal heat abstraction (Fig. 3b). 

The assumption that the process along the path 5-1 is 
isothermal (T = idem) can raise some doubts, as this process 
is difficult to execute. Moreover, we cannot determine the 
internal energy in this process, as for T = idem the specific 
thermal capacity c., — oo. But there is no logical, nor empirical 
reason for not using this process as the reference process. Its 
applicability is justified by the facts that: 

* although in time t — t, — t; (when the cylinder is filled up 
with the fresh air) the heat is delivered to the fresh charge, 
which is heated by the elements composing the working 
space, but this time is short and the temperature does 
not change sufficiently to make adopting the assumption 
T, = T, = T = idem impossible. In cases when the 
temperature is subject to small but recordable changes 
we can assume that the heat supply to the fresh charge is 
the isothermal process taking place at a constant average 
temperature T, ,, 

e since the time t of this process is short, we can assume that 
in that time the fresh charge mass (m) and consequently 
its internal energy (U) do not change remarkably, i.e.: if 
m = idem then U = idem at T = idem. 

A distinctive advantage of this process is that it is most 
economical. 

For the modified Diesel cycle taking into account first the 
isochoric and then the isobaric heat abstraction, the theoretical 
work done in the cylinder during one cycle is given by the 
formula [7]: 
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S, - S;- idem 


Fig. 3. Modified theoretical cycles: a) Diesel cycle taking into account 
isobaric heat abstraction, b) Diesel cycle taking into account isothermal 
heat abstraction: L,, — theoretical work of the cycle, p — pressure, 

p, - average theoretical pressure, V — volume, T — temperature, S — entropy, 
Q,- supplied heat, Q, — heat abstracted at constant volume (isochoric), 
Q, — heat abstracted at constant pressure (isobaric), Q,, — heat abstracted 
at constant temperature (isothermal), GMP and DMP — upper and lower 
dead centre of the piston, respectively 


La =Q - Q», +Q, )= 


(8) 
= me, [«(T, -T,)- (T, -T, + «(T,-T,))] 
or: 
La =Q,- Q, =p;V; + 

4 2 (9) 
+ [pav -|p;V; -p, V + [pav +p,V, 

3 1 

Pi = P; 


At the same time for the modified Diesel cycle taking 
into account first the isochoric and then the isothermal heat 
abstraction, the abovementioned theoretical work is given by 
the formula [7]: 


La = Q, =Q *Q4)- 


-me, tr, -7,)- (T, -1)]-5,6, -s). °” 


Or: 
La 2Q,-Q,-2p,V,* 
4 1 2 
+ fpav — [pav t [pav +p,V, 


3 5 1 


(11) 


Obviously, taking into account the number of realised 
cycles n and the number of cylinders k, the theoretical work 
(L,) needed for calculating the theoretical engine operation can 
is given by the formula [7]: 


L-kmnL, (12) 


or, after calculating the average theoretical pressure p, (Fig. 
4a), by the formula [7]: 


L, =k-n- L =k-n-p,-(V;-V,)= 


(13) 
=k-n-p,-AV;; 
where: 
V; x Vi =AV;, ME 
k — number of cylinders in the engine, 
n — number of realised theoretical cycles in the operating 
time interval [0, t], 
L, — work ofa single theoretical cycle in the cylinder of 
an ideal engine, 
P, — average theoretical pressure, 
V — cylinder displacement volume, 


V. V, — cylinder working space volumes (having the 
interpretation resulting from the symbols in Fig. 3). 


Comparing the effective work done by the engine in time 
t and the theoretical work which can be done by an (ideal) 
theoretical engine in the same time we can estimate the degree 
of excellence (£, o) of energy conversion into effective work 
(L.) by the examined engine, which can be named the degree of 
excellence of engine operation. This degree can be calculated 
for two cases. The first case takes place when we cannot assume 
that during the time t of engine operation the same effective 
work is generated in each engine cylinder. Then the degree of 
excellence (£, o) of energy conversion in the engine can be 
given by the formula [7]: 


m n 
Ee o a 
— k=l j=l 14 
€ = Á 
Ln(ob) L, 
where: 
La 7 effective work of the j-th cycle in the k-th engine 
cylinder, 
L, — theoretical work calculated using formula (13) for 


the modified Diesel cycle. 


The degree of excellence of energy conversion into 
work, defined by formula (14), describes haw far the energy 
conversion into work in the real engine differs from the 
theoretical work of an ideal engine, after n cycles. 

The second case takes place when we can assume that in 
time t of engine operation the same effective work is generated 
in each cylinder. Then the abovementioned degree of excellence 


(Einon) Can be calculated using formula [7]: 
D. d 
ET n(ob == > L, =knL,, 

d L, bi vod 
where: 
k — number of engine cylinders, 
n — number of realised theoretical cycles in the operating 

time interval [0, t], 

L, — effective work of a single real cycle in the engine 


cylinder. 


Of some interest is how far the operation of the real engine 
as a whole differs from the operation of an ideal engine. 
Analysing the same cases as for the engine degree of excellence 


(Einon) Of energy conversion into effective work (L), we can 


define the degree of engine operation excellence (£p). In the 
first case when we have to assume that different effective work 
is generated in each cylinder in particular time intervals t... 
composing the overall time t of engine operation, the degree of 
engine operation excellence (e, , ,.) can be calculated as [7]: 


Dn(ob) 
m n 
» > Loto 
k=l j=l 
£ = (16) 
Dn(ob) Li 
where: 
ak - effective work of the j-th cycle in the k-th engine 
cylinder, 
tx T time of generation of work Eije 
L, — theoretical work calculated using formula (13) for the 


modified Diesel cycle. 


In the second case when we can assume that in the time t 
of engine operation the same effective work L , is generated in 
each cylinder, the excellence degree (€p) can be expressed 
using the same relation as the degree of excellence (Erno) of 
energy conversion into effective work (L,) by the examined 
engine (15). Consequently: 

n 


2 tox =l 
i=l 

Therefore we can assume that in the second case the 
excellence degree (€p) of operation of an internal combustion 
engine is equivalent to the excellence degree (e...) of energy 
conversion into effective work (L,) in this engine. 

If we take into account the indicated work L, in these 
analyses, we can determine the coefficient (degree) of engine 
operation dissipation, caused by the energy lost for overcoming 
mechanical resistances. The abovementioned degree of engine 
operation dissipation could be defined as [7]: 


m n 
2, 2, Lenton 
Spn(ob) = m3 — — — (17) 
2 Liortox 
k-l j=l 

If we can assume that in the time t of engine operation 
the case takes place that: L, = idem and L, = idem, then the 
engine operation dissipation degree is equal to its mechanical 
efficiency. 

Formulas (9) and (11) can be only used when the paths 
of thermodynamic processes of the cycles are known. On the 
other hand, when we know the functional relations between the 
pressures and the volumes which are characteristic for these 
processes, the above formulas take the forms: 


La = Q, =O, =p3V; + 
V, V, (18) 
+ [p(V)dV -| p, V; - p, V, [pC dV + pV, 


Vs V 
L,-2Q,-Q,-p,V, + 
V, V V (19) 
+ [p(V)dV-| |p(V)dV + [p(V)dV » p; V, 
V3 Vs M 


The above analysis shows that in the case when the 
functional relations p — f(V) for the examined processes are 
known, formula (18) can be used instead of formula (9) and 
(19) instead of (11). 
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REMARKS AND CONCLUSIONS 


— The operation of the internal combustion engine was 
interpreted as the delivery of the required energy in given 
time, which can be expressed in the form of a physical 
quantity having the measure unit called joule-second. 

— When discussing energy related advantages of internal 
combustion engines, we should analyse their operation and 
not only their work. Beside the work itself, the operation 
analysis takes into account also time of its realisation. 

— The proposed method can be applied for calculating the 
theoretical standard operation, which enables comparison 
between different real diesel engines revealing different 
levels of wear. The theoretical work needed for calculating 
this operation was calculated based on the modified 
theoretical Diesel cycles, different for different engine 
rotational speeds. The modification introduced in those 
cycles consisted in additional taking into account the 
heat abstraction according to the isobaric or isothermal 
processes. For low-speed engines the additional heat 
abstraction modelled as the isobaric process was proposed, 
while that in the medium-speed engines was assumed to 
take place according to the isothermal process. 

— The presented analysis reveals that it is not advisable to 
use traditional and well-known theoretical Diesel cycles 
for determining the theoretical work making the basis for 
calculating the engine operation value. 

— The article presents the deterministic method applied for 
evaluating the theoretical operation of the piston internal 
combustion engine used as the main engine on a ship, but 
this method can be applied to any piston internal combustion 
engines. 

— Aseparate problem is developing a method for assessing the 
operation of internal combustion engines in the stochastic 
formulation. In this case we should take into account that 
the energy conversion processes taking place during engine 
operation have stochastic nature. 

— Asimilaranalysis can be performed for four-stroke internal 
combustion engines, and similar formulas can be derived 
for determining their theoretical operation. 

— Thearticle also illustrates the applicability of the assessment 
of theoretical internal combustion engine operation, done 
in the above way, for evaluating real operations of this type 
of internal combustion engines. 
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Exhaust gas temperature measurements 
in diagnostic examination of naval 
gas turbine engines 


Part II 
Unsteady processes 


Zbigniew Korczewski, Assoc. Prof. 
Gdansk University of Technology 


ABSTRACT 


The second part of the article presents the results of operating diagnostic tests of a two- and 
three-shaft engine with a separate power turbine during the start-up and acceleration of 
the rotor units. Attention was paid to key importance of the correctness of operation of the 
automatic engine load control system, the input for which, among other signals, is the rate 
of increase of the exhaust gas flow temperature. The article presents sample damages of the 
engine flow section which resulted from disturbed functioning of this system. The unsteady 
operation of the compressor during engine acceleration was the source of excessive increase 


of the exhaust gas temperature behind the combustion chamber and partial burning of the 


turbine blade tips. 


Key words: technical diagnostics; naval turbine engines; exhaust gas temperature; 
unsteady processes; unstable compressor operation; partial blade tip burning 


INTRODUCTION 


In unsteady states of turbine engine operation, when it 
changes from one steady-load state to the other, operating 
parameters and power indices ofthe thermal and flow processes 
taking place inside the engine, as well as the processes 
themselves change in time. This happens, for instance, when the 
engine is started up or shut down, during engine acceleration!) 
or deceleration?, as well as during changes of direction of 
rotation of the driving turbine rotor — in reversal engines. 

The behaviour of the naval turbine engines in unsteady 
states makes their users perform specific actions to secure 
long-lasting and faultless engine operation. On the one hand, 
the user should know what course the transient processes 
take, how the operating parameters change, where and when 
the overload phenomena take place and what is their nature 
and time of duration. On the other hand, the opportunities for 
performing necessary experimental tests of engines in marine 
conditions are strictly limited by numerous structural and 
operational restrictions. 

The operating practice of naval engines reveals that the 
realisation of transient processes in extremely unfavourable 
conditions may result in exceeding the exhaust gas temperature 
limits and unstable operation of the compressor and the 


combustion chamber, which may finally lead to partial burning 
of turbine blade tips [1, 2, 6, 9]. This refers in particular to the 
engine start-up and acceleration processes. 


ENGINE START-UP AND ACCELERATION 


The start-up of the naval turbine engine 1s an unsteady 
process oriented on bringing the rotor units from the state of 
rest to the idling rotational speed. It is mainly characterised by 
the fact that the largest changes of energy conversion processes 
are recorded during its realisation. In as little as a few seconds 
the structural components of the flow section are subject to 
strong thermal loads, with the instantaneous temperature of 
the inflowing exhaust gas increasing by AT = 60 — 70 K/s, 
Fig. 1. This excessive increase of the exhaust gas temperature 
may be provoked by the worsened technical state of the flow 
section, or result from faulty operation ofthe automatic engine 
load control system which controls the rates of fuel and supply 
air [4, 5, 7, 11]. 

Fig. 2 shows sample time-histories of the exhaust gas 
temperature behind the low pressure turbine and behind the 
exhaust gas generator as the functions of the start-up process 
time duration. These time-histories were recorded during 
five years of operation of a three-shift engine with a separate 


D acceleration - process of increasing the rotational speed of the rotor units. 
P? deceleration - process of decreasing the rotational speed of the rotor units. 
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Fig. 1. Time-histories of parameters controlling the start-up process of the three-shaft engine with a separate power turbine nSNC, nSWC — rotational speed 
of, respectively, the low pressure and high pressure compressor; Ir, Ur, Uw — current parameters of the electric starter: strength and voltage of the starting 
current, and the field voltage; T4 — averaged temperature of the exhaust gas flow behind the exhaust gas generator; p2 — air pressure behind the high 
pressure compressor, ppal — fuel pressure in front of the injectors 
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Fig. 2. Exhaust gas temperature changes recorded behind TNC of a three- 
shaft engine with a separate power turbine vs. start-up time: 1 — brand-new 
engine; 2 — engine in third year of operation; 3 — engine after five years of 
operation; 4 — unsuccessful engine start-up after five years of operation; A 
— automatic engine shut-down 


power turbine on a ship. A characteristic trend is observed for 
continuous increase of the exhaust gas temperature derivative 
when the engine's operational time increases. 

The curves shown in the figure reveal that the increased 
worsening of the technical state of the engine flow section 
(confirmed by systematic endoscopic examinations [9]) has 
finally led to the increase of the exhaust gas temperature to 
the limit at which the automatic control system stopped the 
start-up process, thus protecting the engine against the damage 
(point A on curve 4). 

An additional threat for the reliability of the started-up 
engine is possible excessive deflection ofthe cooperation line on 
the compressor characteristic from the optimal line and entering 
the zone of unstable operation (so-called compressor stall) 
[3, 8, 12]. The graphical interpretation of the both undesired 
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> 
m 
Fig. 3. Changes of range of cooperation between compressor, combustion 


chamber, and turbine during engine start-up: 1 — cooperation line 
in steady-load states; 2 — change of operating range at dT’ /dt = dop; 
2'— change of operating range at dT" Jdt = kryt; 3 — stable operation limit; 
ny? ng, — isodromes of idling and electric rotation of the compressor- 
turbine unit, WW, — lines of constant fuel charge for idling 


and electric rotation (beginning of combustion); T',. » T ,,,— isotherms 


of the exhaust gas flow leaving the combustion chamber 


phenomena is given in Fig. 3 which shows the combined 
characteristics of the compressor-combustion chamber- 
turbine system, with projected hypothetical curves of system 
cooperation during the engine start-up process. 

The fuel ignition and the increasing exhaust gas temperature 
shift the points of system cooperation from the optimal line 
for steady-state conditions towards the stable operation limit, 
which is illustrated by lines 2 and 2’ in Fig. 3. The margin of 
stable operation of the compressor, defined as: 


AZ, = (Z, - 1)100% (1) 
where: . 
— (Ts /m),. Q) 
° (Ts /m), 
Z. — stable compressor operation margin coefficient; 


(ng/ m), — ratio ofthe compressor's compression to the mass 
flow rate of the working medium at the stable 
operation limit; 

(n s/m), — ratio ofthe compressor's compression to the mass 
flow rate ofthe working medium when the engine 
switches from the stable speed of electric rotation 
to the stable speed of idling, is reduced. 


Reducing the stable operation margin ofthe turbine engine 
compressor below 5% results in boundary layer separation on 
convex surfaces ofthe rotor blades and pulsations ofthe forced 
air flow [3, 8, 10]. 

Similar gasodynamic phenomena take place in the process 
of acceleration of the engine rotor units when the engine 
shifts from one stable-load state to the other. It results from 
the equation of motion of each rotor unit that the acceleration 
time t, depends on the mass polar moment of inertia J of its 
rotating system and the range of rotational speed changes from 
n, to n,, according to formula [3]: 


CO " 


tT =— 
" 30 -M 


ni p 

Increasing the values of the accelerating torque M, of the 
engine shaft, and/or reducing the values of the mass polar 
moment of inertia of the rotating masses and the rotational 
speed range shorten the time needed for changing the operating 
engine parameters to the required higher-load regime. Therefore 
to increase the engine power, the rate of the fuel supplied to the 
engine combustion chamber for each intermediate rotational speed 
ofthe rotor unit should be larger than in case of engine operation 
in steady-load conditions at the same rotational speed. 

The extremely intensive course of the engine fuel supply 
process during the acceleration Mal = f(n) is limited by the 
maximal permissible exhaust gas temperature measured at 
high-pressure turbine inlet, due to the creep-resistance of 
constructional materials of the turbine blades, the efficiency of 
operation of the turbine cooling system, and the limit of stable 
operation of the compressor. The maximal permissible excess fuel 
charge determines the minimal time of engine acceleration. 

The reserve of stable gasodynamic operation of the 
compressor AZ, which is used in the acceleration process is 
calculated for each instantaneous engine rotational speed as the 
product of two coefficients: AZ — the reserve of compressor 
stability determined by the exhaust gas temperature at turbine 
inlet when the engine operates in steady-load conditions, and 
AZT the reserve of compressor stability which characterises 
permissible exceeding of the exhaust gas temperature at turbine 
inlet during acceleration: 


* 


Tq 


= _ spala spalmax __ 
AZ, = AZ AZ, = T ` T (4) 
spal max spalu 
where: 
T — exhaust gas temperature at turbine inlet in the 


spalmax 


maximal steady-load regime, 

— maximal exhaust gas temperature at turbine inlet 
during the acceleration process, 

— exhaust gas temperature at turbine inlet in the steady- 
load regime. 


* 


spala 


* 


spalu 


In order to protect the engine against damages in these 
situations, the procedure of its use includes numerous limits 
which are automatically controlled by the start-up and load 
change process control system. The role of this system is to act 
on the rate of the energy delivered to the engine (as a result of 


burning the fuel supplied from the fuel installation and from the 
electric starter) and control the flow duct geometry, including 
the adjustable inlet stator, air bleeding valves behind SWC, 
supply air valves for the second injector channel, reversal 
mechanism of the separate power turbine, etc., in order to 
protect the engine against excessive thermal and mechanical 
overloads of its constructional components. In Fig 3 this range 
is represented by the hatched area between lines 1 and 2. The 


block diagram of the system is shown in Fig. 4. 
Fig. 4. Block diagram of the turbine engine exhaust gas temperature limit 
control system: A — impulse temperature control module, B — module 
of actuating mechanisms, C — naval turbine engine (OTSS), T — module 
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of thermo-couples (TXA-1368), Z — temperature setter, R — computer 
programme for measuring, recording and analysing engine start-up 
parameters, 1 — signal amplifier and modulator, 2 — electromagnetic 
overflow valve (by-pass) control circuit, 3 — main valve and emergency fuel 
drainage valve control circuit, 4 — main fuel valve (ZGP), 5 — emergency 
fuel valve (ZAP), 6 — fuel overflow valve (ZPP), 7 — comparing element, 
8 — reference voltage generator, e — averaged thermoelectric voltage value, 
s — pulse control signal 


The input signal for the system is the average exhaust 
gas temperature measured behind the exhaust gas generator. 
Excessive increase of this temperature can “only” provoke the 
action of the fuel overflow valve, or, in more severe situations, 
rapidly shut the engine down by acting simultaneously on 
the main valve and the emergency fuel drainage valve. 
It is noteworthy that the measurement of the exhaust gas 
temperature in the start-up process (and other unsteady 
processes) is burdened with some relatively large inertia 
which can be corrected electronically in the engine control 
system — Fig. 5 [13, 14]. Correction of the measuring inertia 


T 


Fig. 5. Hypothetical time-histories of exhaust gas temperature changes 
during turbine engine start-up (acceleration) — the hatched area represents 
the range of correction of the measuring thermocouple inertia: 

1 — real temperature, 2 — measured temperature 
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Fig. 6. Partly burned high-pressure turbine rotor blade tips 
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Fig. 7. Time-histories of high-pressure rotor rotational speed, negative pressure at compressor inlet, and exhaust gas temperatures 
at high-pressure turbine and separate power turbine inlets during the two-shaft engine acceleration process 
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Fig. 8. Changes of exhaust gas temperature at high-pressure turbine inlet during start-up, acceleration, 
and steady-load engine operation of the two-shaft engine with a separate power turbine vs. high-pressure rotor rotational speed 


of the thermocouples, which is done individually for each 
engine due to unrepeatability of the process of engine 
manufacturing, may reach as much as a few percent. This 
correction determines a so-called starting thermal sensitivity 
of the engine. On the other hand, interfering into producer's 
settings of the automatic control system during all time 
of engine operational use 1s not allowed, in particular in 
situations when the engine start-up becomes relatively 
difficult due to the worsened technical state of the flow section 
or the fuel supply system [13]. 

The absence or disturbed functioning of the automatic 
exhaust gas temperature limiter usually leads to huge and 
irreversible damages ofthe flow section, after which recovering 
the technical state ofthe engine is only possible via emergency 
repair done by the producer. Fig. 6 shows partly burned blade 
tips of the high-pressure turbine rotor, being the result of 
unstable operation of the combustion chamber during the 
acceleration process of a two-shaft engine with a separate 
power turbine. 

After accidental choking of the air inlet duct and closing 
of the air bleeder valve, the increased rotational speed of 
the exhaust gas generator rotor was accompanied by rapid 
increase of the negative pressure in the compressor's inlet 
section, the decrease ofthe compressor output and, as a further 
consequence, by the increase of the exhaust gas temperature 
behind the combustion chamber — Fig. 7. The maximal 
permissible value of the exhaust gas temperature in steady- 
load states, limited to 970°C, was considerably exceeded, as 
shown in Figs. 7b,c and 8 (exceeding of the upper limit of the 
measuring range). Another parameter which was considerably 
exceeded in this case was the exhaust gas temperature at the 
separate power turbine — Fig. 7d. 

In unsteady processes some short-lived exceeding of the 
maximum exhaust gas temperature at high-pressure turbine 
inlet is permissible, but this temperature is definitely not 
allowed to increase by over 2096 of its maximal value referring 
to the nominal steady-load engine operation (depending on the 
used constructional materials and applied method of turbine 
cooling) [3, 6]. 

The excess rate ofthe fuel burned in the combustion chamber 
has led to excessive enrichment of the fuel-air mixture and the 
decrease of the excess air number X. — Fig. 9. When the critical 
level was reached (A << Amin) ata given rotational speed of the 
exhaust gas generator rotor unit, the phenomenon of so-called 
rich" flameout of the combustion chamber (flame breakdown) 
took place, with further automatic shutting the engine down. 
But despite this shut-down, the temperature ofthe flow section 
components (high pressure turbine blades) was so high that the 
resultant damages (partial burning, cracks) ofthe constructional 
material have eliminated the engine from further operation. 
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Fig. 9. Relation between the excess air number in the combustion chamber 
and the rotational speed of the exhaust gas generator rotor unit 


CONCLUSIONS 


— The rate of increase of the exhaust gas temperature during 
the turbine engine start-up and acceleration processes is 
decisive about engine reliability and lifetime. It is also a 
basic diagnostic parameter monitored in unsteady processes 
in order to assess the technical state of the flow section, the 
fuel and air supply system, and the automatic load control 
system. 

— Akey operational problem is adjusting the engine in a way 
which will simultaneously secure sufficient dynamics of the 
realised unsteady thermal and flow processes and, on the 
other hand, protect the engine against exceeding the limits 
of the exhaust gas temperature and the rate of its increase 
as possible threats of vast and irreversible damages of the 
engine flow section. 

— Therefore one of basic requirements for the thermal 
protection system of the engine is maximal reliability of 
its work during engine operation, which is checked in 
systematic diagnostic tests. 
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Damping of torsional vibrations of ship engine 
crankshafts — general selection methods 
of viscous vibration damper 
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ABSTRACT 


This paper describes causes of torsional vibrations generated in ship engine crankshafts. 

Means for damping the torsional vibrations as well as general methods for selection 

of viscous torsional vibration dampers for a given type of engine are also presented. 

Exemplary results of calculations connected with selection ofa viscous torsional vibration 
damper intended for a six-cylinder engine, are attached. 


Keywords: torsional vibrations; torsional vibration damping; torsional vibration dampers; 
crankshaft; natural vibration frequency 


INTRODUCTION 


Crankshafts are applied everywhere it is necessary to 
convert reciprocating motion into rotational one or inversely. 
In engines and compressors they are commonly used. 

Piston combustion engine in operation generates vibrations 
which result from occurrence of periodically varying gas and 
inertia forces. The forces generate the following kinds of 
vibrations: 

* bending vibrations, 

* axial vibrations, 

* torsional vibrations. 


The first two kinds of vibrations occurring in car combustion 
engines do not constitute any great danger for life time of their 
crankshafts but for ship large-power engines axial vibrations of 
crankshafts are a serious problem. Such vibrations cause that 
the entire system consisted of engine crankshaft, flywheel, shaft 
line and screw propeller displaces periodically along its axis. 
Axial vibration amplitude of the system practically depends 
on a design solution of screw propeller, namely, number of 
its blades. 

However irrespective of an used dynamical engine 
operation system torsional vibrations are the most dangerous 
for its crankshaft. Out of many forces acting in the crankshaft- 
piston system the force T (Fig. 1) tangent to the circle described 
by crank, makes crankshaft rotating. The force T, one of two 
components of the force S acting along crankshaft axis, is that 
periodically changing [8, 9, 11, 12]. Its frequency for two-stroke 
engines is equal to 2x, and for four-stroke ones - to 4z. Changes 


Fig. I. Distribution of forces in crankshaft - piston system 
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of the force T generate accelerations of engine crankshaft 
rotational motion and in consequence torsional vibrations which 
change along with changing shaft line rotational speed. Worth 
mentioning that crankshaft deflection is limited only by shaft 
line torsional stiffness and torsional vibration amplitude can 
exceed its permissible values. In the case of lack of damping 
vibration run curve tends to infinity for any speed which 
generates a resonance harmonic. 

Torsional damage of crankshaft, i.e. one of a changeable 
stiffness, takes place when the permissible amplitude value 
(permissible twist angle) o ie is exceeded (1) [2, 3]. 


9. > P dop 
2L. (1) 
Qa = 

op s = Gly 
where: 
M, — torque, 
L, — reduced length of shaft segment, 
G — shear modulus, 
I. — polar moment of inertia. 


oi 


Many formulae are used for calculating the reduced length 
of crankshaft, i.e. shaft segment between cylinder axes. Among 
the most known the following may be numbered: 

e Carter formula, 
* Ker Wilson formula, 
* Tulpin formula. 


The formulae may be reduced to the general form as follows 
(2) [1, 2, 3, 8, 9]: 


L,=k, lk, leak, 22 (2) 
i cg I ck I TW I 
where: c S id 
I — cross-section inertia moments, 
k — set of quantities which account for crank geometry, 
and the indices stand for: 
cg — main journal, 
ck — crank journal, 
rw-— crank arm, 
o — reduced value. 


Importance of the problem in question, 1.e. danger of shaft 
damage due to occurrence of torsional vibrations in propulsion 
systems of floating objects resulted in that Polish Register of 
Shipping (PRS) issued a specialty rule publication entirely 
devoted to the calculating of crankshafts of combustion- 
ignition engines [10]. One of its sections contains guidelines 
and conditions to be abide by crankshaft designers. 


DAMPING OF TORSIONAL VIBRATIONS 
OF CRANKSHAFTS OF SHIP ENGINES 


For many years for the damping of torsional vibrations of 
ship engine crankshafts the following kinds of dampers have 
been used: 

e frictional dampers, 

* rubber dampers (Photo 1), 
* viscous dampers (Photo 2), 
* spring dampers (Photo 3). 


Among the above mentioned group of dampers the viscous 
one is most popular and successfully used for controlling 
vibration level of ship large-power medium-speed engines 
[4, 5]. 
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Photo 1 


Photo 2 Photo 3 


GENERAL METHODS OF SELECTION AND 
DESIGN OF VISCOUS DAMPERS 
OF TORSIONAL VIBRATIONS 


In advance of commencing design of a viscous torsional 
vibration damper its designer should obtain, from engine 
producer, an appropriate set of necessary data [6]. 

The data should contain the information on: 

* kind of fuel (petrol, diesel oil), 

* number of cylinders, 

* type of engine (two-stroke, four-stroke), 

* minimum rotational speed of engine, 

* maximum rotational speed of engine, 

* range of operational speed of engine (if not constant), 
* ignition sequence, 

e main journal diameter, 

* crank journal diameter, 

* cylinder diameter, 

* mean indicated pressure, 

* total weight of masses in reciprocating axial motion, 
* torsional stiffness of shaft between cylinders, 

* permissible value of twist angle of shaft, 

* limitations for damper gabarites. 


The above specified data make it possible to form, for real 
propulsion system, a substitute vibrating system which enables 
to calculate: 

- its natural frequency Q [rad/s] (3), 
- resonance rotational speeds: 


| 600. 


D. (3) 
wh 2m h 
where: 
Q — natural vibration frequency of w-form, 
w — vibration form number, 
h -  harmonic's order; for two-stroke engines: 


h=1, 2,3, 4,...; 
for four-stroke engines: 
h=0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4... 


- relative torsional vibration amplitudes ©, of particular 
masses of the system (Fig. 2), as well as: 

- range of harmonics in the area where calculations have to 
be performed (4): 


E 609, 
2nn 


max 
min 


60Q.. 4) 


min 
200... 

For the so determined range of harmonics, values of the 
tangent force T, are calculated with accounting for the mean 
indicated pressure. Despite that the determined range of 
harmonics may be very wide it is sufficient, basing on multi- 
year experience, to take into account, in the performed analysis, 


only a certain number of harmonics of the tangent force T. It 
is usually 12 +18 harmonics, because higher-order harmonics 
of wide ranges and small amplitudes are of low impact on 
vibrations. 
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Fig. 2. Exemplary diagram of relative amplitudes calculated for a six- 
cylinder engine; 1) damper, 2 + 7) crankshaft- piston system, 8) flywheel 


As the following is known: 

* values of relative amplitudes of masses, 

* range of harmonics, 

* values of tangent forces, 

* values of rotational speed, which correspond to a given 
harmonic, 


torsional vibrations of crankshaft free end can be calculated. 
First of all the calculations are performed without accounting 
for torsional vibration damper in the system, but with taking 
into account screw propeller and elastic coupling (1f applied). 
Value oftorsional vibration amplitude of engine crankshaft free 
end is described by the following formula (5): 

c 


90A $9, 
©, 2——À—— (5) 
2B+2C+2D 
where: 
ASI KT. 
Bz8,5 070... 
i=l 
E 2 
C=8, 05 Qv 
2 
D = k, 9. y 
M 

©, — torsional vibration amplitude of crankshaft free end, 
D  — cylinder diameter, 
R — crank radius, 
T, — tangent component of exciting harmonic, 
©, — relative amplitude of i-th mass, 
XO, — geometrical sum of relative amplitude vectors, 
c  — number of cylinders, 
S. - engine structural damping coefficient, 
S, — screw propeller damping coefficient, 
k — coupling torsional rigidity, 
M - coupling dynamic amplification (specified by coupling 


producer). 


In order to obtain full image of crankshaft loading resulting 
from torsional vibrations, value of torsional stresses is 


calculated, which can be analyzed by using diagrams prepared 
on the basis of the performed calculations. 

The determination of values of twist angles of crankshaft 
free end and torsional stresses closes the first phase of damper 
selection after which preliminary calculations of the mass 
moment of inertia of damper inertia ring, I , are performed. 

The mass moment of inertia of damper inertia ring, I , for 
a selected (designed) damper, has to comply with the following 
condition (6): 

I, 


Ye I; 
i=l 


0.4< «I 


(6) 
where: 
> ©? I, — the so called equivalent inertia of engine. 
i=l 
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Fig. 3. Diagram of changes of torsional stresses in function of shaft 
cross-section, for a selected harmonic; exemplary results 
of calculations for a 6-cylinder engine 
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Fig. 4. Diagram of changes of torsional stresses in function of harmonic, 
for a selected vulnerable cross-section of shaft; exemplary 
results of calculations for a 6-cylinder engine 


From the above given equation unambiguously results 
that for further analysis it is possible to obtain practically a set 
of ,,1” dampers whose effective mass moments of inertia are 
expressed as follows (7): 


I, 
I, =I, += (7) 
i i 2 


where: 


I; — mass moment of inertia of i-th damper casing. 
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Now, having the preliminarily selected set of,,i” dampers, one 
performs identical calculations as for the system without damper 
but extended by the effective mass moment of inertia I 

For every i-th system equipped with i-th damper the 
following is determined: 

* system natural frequency, 

* range of harmonics, 

* resonance speeds, 

* values of relative amplitudes (Fig. 5), 

* values of twist angles of crankshaft free end, 
* values of torsional stresses (Fig. 6). 


Comparative analysis of the obtained results makes it 
possible to select, from the set of ,,i” dampers, that most 
effective in a given dynamic system. 
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Fig. 5. Diagrams of relative amplitudes of particular equivalent systems 
for the whole set of preliminarily selected dampers; exemplary results 
of calculations for a 6-cylinder engine 
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Fig. 6. Diagrams of torsional stresses in function of the most loaded 
cross-sections of shaft for the whole set of preliminarily selected dampers; 
exemplary results of calculations for a 6-cylinder engine 


Viscous torsional vibration damper is that to which the fluid 
friction phenomenon is applied to effect vibration damping [7]. 
Therefore also value of friction moment resulting from friction 
forces acting on working surfaces of damper ring and casing, 
decides on its effectiveness. 

The total value of friction moment which appears on damper 
inertia ring may be expressed as follows (8): 
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M, = M + M ei T M,» (8) 


and, that of friction moment on working surfaces of damper 
casing (9): 

M. M oz + M ow + M o (9) 
where: 
M, - 
and the indices stand for: 
p - ring, 

o — casing, 

z — outer surface, 
w 

b 


friction moment, 


— inner surface, 
— side surface. 


Value of the friction moment does not depend only and 
solely on the geometrical parameters of damper but also on 
the fluid kinematic viscosity (10): 


n 


v=— 
p 


(10) 
where: 
n mE 
p = 


dynamic fluid viscosity, 
fluid density. 


During selection of fluid viscosity intended for damper it 
should be remembered that the viscosity changes along with 
temperature changing. 

As already mentioned, the vibration damping in the designed 
device is effected due to friction, hence energy dissipated by 
damper is transformed into heat. The heat generated in damper 
during the initial phase of its work is transferred mainly to 
the environment and partly used for heating the damper. It so 
long proceeds until the damper obtains the so-called saturation 
temperature. 

After reaching the saturation temperature by the damper 
a justified concern arises that it may become seized due to 
its overheating resulting not only from the transformation of 
dissipated energy into heat but also its operation in an elevated 
environmental temperature (operation of damper in engine 
chamber) (Photo 4). 


Photo 4. Image of an overheated and seized viscous torsional vibration 
damper, acc. tests made by the DAMPOL 


In practice the detrimental phenomenon makes that the 
damper, instead to damp vibrations, starts to excite them. 
The damper, if correctly designed, should operate reliably 
for 20000 + 30000 h in negative temperatures and first of all 
higher temperatures even as high as 120°C. Laboratory tests 
performed on viscous torsional vibration dampers demonstrated 
that events of seizing the dampers occurred when their working 
temperature values exceeded by 60°C the temperature for which 
they were designed. 

Satisfactory effects of damper operation are reached when 
it operates within the temperature range from 75° + 90°C in 
continuous cycle of work. 


In view of damper operation effectiveness it is recommended 
to keep the specific power dissipated by damper lower than 
8.6 [KM/m?] (counted per total surface area of its inertia 
member). 

The range of heat flow within damper depends on: 

* working conditions (a.o. engine rotational speed) 

e damper size, 

* damper design (of screw bolted casing or cold shut one) 
and is expressed as follows (11): 


Heh-A. (11) 


where: 

H - permissible heat flow rate [J/h], 

A, — total surface area of inertia member, 

h = 18 + 22 [MJ/m?/h] — for instantaneous operation at critical 
speed, 

h=9 + 11 [MJ/m?/h] — for dampers used for small high-speed 
engines intended for continuous operation at critical speeds, 
such as those applied to cars, 

h = 4.5 + 5.5 [MJ/m?/h] — for dampers of large low-speed 
engines intended for continuous operation at critical speeds, 
such as those used in ships. 


Note: The paper has been financially supported by the 
National Science Centre - Research Project No N N509 
547440. 
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ABSTRACT 


The article presents the operational and economic analysis of controlling the cooling water flow in marine 

steam turbine power plants. The analysis bases on selected designs of the main condenser cooling water 

pumps and makes use of the results of investigations performed in inland power plants. Special attention 
was focused on marine aspects of the operation of those systems. 


Keywords: steam turbines; turbine operation; marine power plants; marine pumps; circulating pumps 


INTRODUCTION 


One of basic components of the steam turbine power plant 
is the condenser in which the steam condenses after passing 
its heat to the cooling water. In marine power plants the water 
used for cooling the steam is taken from the outside, i.e. from 
the water area in which the ship sails. In the majority of power 
plants the condenser cooling water pumps work at constant 
mass flow rate, but recently a more and more frequently 
analysed problem is an attempt to control this parameter, 
especially when the turbine works under changing load 
conditions. The main reason for applying control systems in 
inland turbine power plants is an attempt to save the water, 
the consumption of which is extremely high during steam 
power plant operation [1, 5, 8, 12]. Marine turbine power 
plants also reveal remarkable load changes. Although turbine 
driven tankers for transporting liquid oil-based chemicals 
have not been built any longer in recent years, a number of 
such tankers for transporting liquefied gas were launched. 
Therefore the economic analysis of cooling pump operation 
in marine power plants is fully justified, especially in case of 
tankers, which sail with the cargo in forward trip and return 
only loaded with the ballast. The above situation considerably 
differentiates ranges of power required by the ship propulsion 
system. Moreover, ships with unlimited sailing region, which 
are, as a rule, designed for least favourable conditions of 
power plant operation in the tropical zone, also sail in more 
moderate zones, which needs adapting the operating conditions 
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of their cooling systems to current ambient temperatures. 
Another argument justifying the interest in analysing the 
above problem is the appearance of opinions suggesting the 
need for reconstruction of the global fleet of liquefied fuel 
tankers based on new standards [17] due to the exploitation 
of new deep-sea oil development areas. We can expect that in 
the nearest future new tankers with both engine and turbine 
propulsion systems will be built. 

Until recently, attempts to decrease the power used by 
the auxiliary devices have been almost unprofitable, which 
was connected with the loss of turbine unit power due to the 
pressure increase in the condenser caused by the reduced mass 
flow rate of the cooling water at lower loads. It was partially 
justified, as the efficiency gain resulting from such an approach, 
which was estimated based on on-line measurements, turned 
out very inaccurate. 

That is why the mass flow rate of the condenser cooling 
water was controlled rather rarely. This control referred to, 
for instance, changing the angular position of the blades in 
the rotors of axial-flow pumps, done in order to reduce the 
technically unfavourable power limits during pump start-ups 
rather than for economic and water saving reasons [15, 16]. 

It should be mentioned here that the “seasonal” rotation 
of axial-flow pump rotor blades, which was used for years, 
resulted from the "numerical" water flow optimisation 
oriented on adapting the system to tropical, moderate, 
or arctic conditions. The above actions were a source of 
operating problems connected with blade seizing in these 
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Fig. 1. Designs of cooling water pumps used in inland 360 MW steam turbine power plants (a) and their rotors used in 200 MW power plants (b) [2, 3] 


pumps. Pumps of other types have no practical possibilities 
to control the water flow via changing the geometry of the 
rotating elements. 

The article presents a comprehensive evaluation of the effect 
of controlling the mass flow rate of the water pumped by the 
diagonal pumps for cooling the condenser on the economy of 
operation of the entire power plant. The analysis takes into 
account, on the one hand, the reduction of the power needed 
for driving the cooling water pumps at their reduced capacity, 
but also, on the other hand, power losses of the main turbines 
resulting from the lower-quality vacuum in the condenser due 
to smaller rate of the cooling water. Operating aspects of such 
solutions which are specific for marine applications are also 
discussed. 


PROPOSED DESIGNS OF COOLING 
WATER PUMPS 


Controlling the cooling water mass flow rate started to be 
the subject of discussion when the turbine blocks in inland 
steam turbine power plants were to be modernised. Old pumps 
in those blocks were subject to failures during block operation, 
as they frequently worked in cavitation areas. Consequently, 
they did not guarantee preserving the required parameters in the 
repaired blocks. Castings of their rotors were made carelessly of 
carbon steel, which considerably affected their lifetime. Modern 
design and technological procedures and methods have made 
it possible to design new rotors which remarkably improved 
the efficiency of pumping and provided opportunities for better 
adjustment of the operation of these pumps to changing load 
conditions. The designs of the cooling water pumps for inland 
large-power turbines are shown in Fig. 1 [2, 3]. 

Tests of the new pumps have revealed remarkable 
advantages of their designs, which can be observed when 
comparing the operating characteristics of the old and new 
pumps cooperating with one turbine, see Figs. 2 and 3, [7]. 
The efficiency gains after modernisation amounted to over 
10%, at simultaneous increase of the elevation head and water 
mass flow rate. 

Positive tests were the motivation for designing water 
control systems for those pumps to adapt their operation to 
changing power loads and environmental conditions. The 
technical analysis performed for the existing configurations 
of the cooling systems has revealed that the best way to 
control the water mass flow rate is the use of guide vanes. 


Introducing guiding blades at the impeller pump inlet turned 
out operationally stable, cheapest, and very efficient when the 
control range was not excessively high [14]. An inlet guide vane 
designed and installed in the cooling water pump working in 
the 360 MW turbine system is shown in Fig. 4 [4]. By changing 
the guide blade setting angle a the cooling water flow could 
be changed with high efficiency by even as much as 30% 
reduction of the nominal mass flow rate, what can be seen on 
the experimentally determined operating characteristics, Figs. 2 
and 3. This method turned out much better than throttling, 
which was used in the past for pumps cooperating with turbines 
and cooling towers [1]. 

It also turned out that the inlet guide van blades can be set 
with "negative" angles with respect to the direction of rotor 
rotation, thus introducing an opposite angular momentum of 
the water flow onto the rotor blades. Such an approach enabled 
increasing the mass flow rate of the condenser cooling water 
over the nominal level for unchanged pump construction. 
Indeed, this effect is occupied by increased power taken by the 
motor used for driving the pump, but it makes it possible to 
decrease the steam condensing pressure when the temperature 
of the cooling water is extremely high. This pump operation 
control is only possible when the (usually electric) motor 
used for driving the pump has excess power reserve and can 
be overloaded. 
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Fig. 2. Generalised flow characteristics of the cooling water pumps 
installed in an inland power plant: new pump (1 ~ 3), 
old pump (4), pipeline (5), [3] 
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Fig. 3. Efficiency characteristics of the old non-controlled pump (6) 
and the modernised pump (3), along with the electric power consumption 
of the cooling water pumps installed in an inland turbine system: 
old pump (1) and new controlled pump (2, 4, 5) [2, 3] 


COOLING WATER PUMP CONTROL 
IN SEA CONDITIONS 


The operation of the cooling water pumps in marine turbine 
power plants was performed using the gas tanker driven by 
the main turbine of an approximate power output of 25000 
kW as a sample case. A scheme of this turbine power plant is 
shown in Fig. 5. 

The main condenser is cooled with the sea water using 
circulating pumps of approximate capacity 8000 m*/h (about 
2250 kg/s) each. The configuration of the condenser cooling 
system is given in Fig. 6, while the scheme of cooling water 
pump connections is shown in Fig. 7. 

Since the main turbines work at changing power, it 1s 
advisable to use cooling water pumps of changing capacity. 
For this purpose the control systems described in the previous 
section which make use ofthe inlet guide vane can be applied. 
These systems have already been positively tested and checked 
in a number of inland power plants. It is noteworthy that the 
cooling water mass flow rate is controlled by changing the 


Fig. 4. View of inlet guide vane blades (a) and blade setting mechanism (b) 


setting of the inlet guide blades, i.e. the immovable element. In 
inland condition, no traces of seizing of the setting mechanism 
of the guide vane blades were recorded during the pump 
operation. The above control method improves operating 
characteristics as it requires less servicing actions, a factor 
which is very important in sea conditions. 


Fig. 5. Scheme of steam turbine power plant of approximate power of 25000 kW in an LNG tanker [17] 
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Fig. 6. Configuration of the main condenser cooling system in the steam turbine power plant 
of an approximate power of 25000 kW in an LNG tanker [17]; a) cross section; b) view. 


A system feature which should be mentioned here is the 
possibility to increase the mass flow rate of the condenser 
cooling water by introducing the negative angular momentum 
at water inlet using a negative setting angle of the inlet vane 
blades. This solution can be profitable for the marine power 
plants operating in extreme tropical conditions. Another way 
of using this control method may consist in assuming lower 
sea water temperatures that those for tropical conditions when 
designing power plants of the ships with unlimited sailing 
region. 


ECONOMIC AND TECHNICAL ANALYSIS 
OF THE USE OF COOLING WATER 
CONTROL 


For the time being, the economic and technical analysis of 
the use of cooling water control systems in steam turbines to 
follow turbine load changes was only possible based on detailed 
design calculations taking into account the turbine geometry, 


see for instance the study by Dzierzgowski [6]. This situation 
originated from the fact that evaluating possible gains of such 
control should take into account the entire thermal cycle of 
the turbine set, with very accurate evaluation of the operation 
of so-called “cold end” apparatuses including: low pressure 
turbine, condenser, pumps, and the cooling water system. It was 
a rather difficult task as the gains resulting from the use of this 
control were very small compared to the turbine power output, 
while the physical description of thermal and flow phenomena 
accompanying this process, necessary for performing the above 
calculations, were extremely complicated. 

This situation was the reason why the advantage of water 
control could not be uniquely evaluated, even based on 
warranty thermal measurements performed in steam turbine 
power plants [8, 12]. In inland conditions this inability was 
checked on 200 MW blocks, where the gains calculated using 
the nomogrammes obtained from the producers and describing 
the operation of these apparatuses [7] were at the level of 
calculating errors. 
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Fig. 7. Connections of main condenser cooling water pumps in the steam turbine power plant of an approximate power of 25000 kW in an LNG tanker [17] 


In this article, the water control effects are proposed to be 
evaluated using the code DIAGAR [9] which was successfully 
used in thermal diagnostics including repair prognoses for 
individual apparatuses in 200 MW turbine systems in inland 
power plants [10, 11]. High accuracy of calculations done using 
this code results from the fact that the balances obtained from 
the measurements are currently verified by the calculations 
making use of detailed description of the geometry of individual 
apparatuses in the thermal cycle. The performed investigations 
have proved that this numerical “calibration” of the results 
of the measurements done by DCS systems in power plants 
provides opportunities for evaluating the operation of individual 
apparatuses at the accuracy level below +30 kW [10]. The above 
methodology has made it possible not only to evaluate precisely 
the efficiency of particular apparatuses, follow the course of 
their degradation, and prepare a reasonable repair plan, but also 
to analyse possible modernisation and improvement variants, 
including the evaluation of advantages of particular cooling 
water control systems. 

Anumerical test of the above described cooling water pump 
control system was performed for a selected inland power 
plant turboset and for a marine power plant shown in Fig. 5. 
The apparatuses in the frame in Fig. 5 represent the “cold end” 
of the turbine and in the examined case were the subject of 
detailed control. 

The initial version of the diagnostic calculations performed 
for the thermal cycles of the inland and marine power plants had 
omitted cooling water system components [12]. For the purpose 
of the present analysis these components were included in the 
calculations, along with relevant data about their operation, 
which made it possible to calculate their energy parameters 
in changed operating conditions. Beside pumps and pipelines, 
the delivered information also referred to low-pressure turbine 
and condenser components which had been modernised in the 
inland turboset. The calculated characteristics of the power 
plant turboset were verified by comparing them with the results 
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of the measurements done on a real block [3, 11]. The obtained 
results are shown in successive figures, with labels facilitating 
their interpretation, to illustrate advantages of the discussed 
control method. 

The pump characteristic shown in Fig. 3 refers to electric 
power and efficiency changes as a function of the water mass 
flow rate. In this case the gains connected with the reduction 
of the mass flow rate result directly from the savings of the 
electric power needed for driving the pumps. This was not 
so clear for old pumps, whose characteristics (1, 6) obtained 
from measurements are given in Fig. 3, for the case of throttle 
governing. The operating characteristics recorded for the control 
performed by changing the blade settings, which are shown in 
this figure, have been recorded in model tests [2, 3]. 

The characteristic of the low-pressure turbine part shown 
in Fig. 8 reveals that the load change, which is here equivalent 
to the change of the volumetric steam flow rate (G, ' V, ) can 
result in efficiency changes equalling even as much as 5%. For 
the given turboset power and the assumed steam mass flow 
rate Gp, the efficiency gains (or losses) result from specific 
volume V, at the controlled mass flow rate G,, of the cooling 
water. The quantity V, is stabilised for new conditions in the 
condenser as a result of pressure changes behind the turbine: 
V, = f(p,). These parameters can be determined directly from 
the condenser characteristics, which for the tested turbine are 
shown in Fig. 9. They are given in the form of p, = f(G.) for 
different mass flow rates and constant inlet temperature of the 
water. 

The same efficiency changes, in percents, were assumed 
for the analysis of the marine power plant. 

Based on the curves shown in the figures, which correspond 
to turbine and condenser operation at nominal and changing 
loads, it is difficult to evaluate precisely gains resulting from 
the cooling water control. It only becomes possible after the 
calculations of the entire thermal cycle, done by the code 
DIAGAR which makes use ofthe above operating characteristics 


of the cold end apparatuses. The results of these calculations, 
which include the evaluation of power gains and the specific 
heat consumption of the entire thermal cycle as a function of 
the load with and without cooling water control, are shown 
in Figs. 10 and 11. They are given in the form of AN; Aq =f 
( N/No) for inlet water temperatures equal to 15°C and 30°C 
which corresponds to winter and summer conditions of turboset 
operation. The diagrams take into account the assumed 15% 
reduction of the water mass flow rate with respect to the nominal 
rate, which was assumed here at the level of 25000 m?/h. 
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Fig. 8. Efficiency of the modernised low part of the inland turbine vs. 
volumetric flow rate of steam for different mass flow rates 
of steam and different condenser pressures 
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Fig. 9. Condenser pressure changes vs. cooling water mass flow rate for 
different mass flow rates of steam at the assumed inlet water temperature 


The presented diagrams reveal that the 15% reduction of the 
water flow in winter (t, = 15°C) brings gains practically within 
the entire operating range. At 50% load these gains exceed 300 
kW (about 0.2% of the turboset power). It is noteworthy that 
the water flow control for low loads eliminates supercooling of 
the condensate in the condenser, see Fig. 10. This is a positive 
500 
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Fig. 10. Increase of turboset power (specific heat consumption) vs. load 
for the cooling water mass flow rate reduced by 15% 
for winter conditions (t. = 15°C) 


feature as the supercooling effect is unfavourable for the 
turboset efficiency. At the same time in the summer, at high inlet 
water temperatures (t, ~ 30°C) the control brings gains only 
when the load decreases below 50%, while the loss at nominal 
load and 15% water flow reduction reaches 400 kW! In this 
case the cooling water flow should be increased rather than 
reduced. Therefore it is advisable to evaluate precisely once 
again the nominal cooling water mass flow rate in the cooling 
pumps for each turboset after its modernisation. 
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Fig. 11. Increase of turboset power (specific heat consumption) vs. load 
for the cooling water mass flow rate reduced by 15% 
for summer conditions (t. = 30°C) 


TECHNICAL ASPECTS OF COOLING 
WATER CONTROL IN MARINE 
POWER PLANTS 


Preliminary calculations performed for the marine power 
plant, see fig. 5, with the aid of the code DIAGAR confirmed 
the opinion that the nature of efficiency changes and power 
gains resulting from the control of the mass flow rate of the 
water cooling the condenser is similar to that recorded in 
inland turbine power plants. However, in case of the marine 
power plants there are more uncertainties concerning the 
temperatures of the seawater taken for cooling purposes. 
Nevertheless we can conclude that the percent changes of 
power and efficiency observed in inland and marine power 
plants are close to each other. 

That means that in the moderate and tropical zones the 
characteristics realised by the marine power plants will 
be similar to those of the inland power plans working in 
comparable conditions. In particular, in the moderate zone at a 
reduced power plant power we can expect the maximal power 
gains within the limits up to 0.2% of the nominal power for 
the entire thermal cycle, taking into account both the reduced 
power needs for driving the pumps, and the power loss 
resulting from the lower-quality vacuum in the condenser. For 
the marine steam turbine power plant with the main turbines 
of the approximate nominal power of 25000 kW it gives 
about 50 kW of the reduction of the ship’s power plant load. 
Assuming that the specific fuel consumption in this power 
plant is approximately equal to 200 g/kWh, we can expect 
the reduction of the fuel consumption by about 240 kg per 24 
hours. The above assessment refers to the moderate zone. The 
expected gains will be much higher when sailing close to the 
arctic zone. 

In extreme conditions in the tropical zone we should 
increase the mass flow rate of the condenser cooling water to 
reach the deepest possible vacuum. We can easily obtain this 
effect by generating a negative angular momentum via turning 
the guide vane blades. 

In marine conditions axial pumps can be used, as the 
elevation head, measured in mwc, needed for cooling the 
condenser is relatively small. These pumps can also be 
controlled using inlet guide vanes. The power and efficiency 
gains expected in this case can be even larger that for diagonal 


pumps [2]. 
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CONCLUSIONS 


The presented analysis has proved that it is advisable to 
control the cooling water mass flow rate in steam turbine 
power plants to follow turbine load changes. Possibility 
of this control is provided in the designs of cooling water 
pumps installed in recently modernised large-power blocks 
[2, 3]. These pumps reveal cavitation-free operation and 
higher efficiency. Water flow changes can be realised in 
these pumps by changing the angular settings of the inlet 
guide vane blades, to response to accurately measured 
changes of the electric power of the turboset and the inlet 
temperature of the cooling water, or (alternatively) the 
pressure in the condenser. However, these signals are to 
be correlated with the calculations of the turboset thermal 
cycle, which can be done online using the data recorded by 
the power plant measuring system and a relevant computer 
code [13]. 

Preliminary calculations making use of the results of 
measurements performed in inland power plants have 
revealed that in winter for the 200 MW blocks, power gains 
can be obtained within the entire range of turbine operation 
when the mass flow rate of the cooling water is reduced by 
15% with respect to the nominal value. When the load is 
reduced by 50%, these gains can even reach 400 kW i.e. 
0.2% of turboset power, which is equivalent to the decrease 
of the specific heat consumption by 15 kJ/kWh. In the 
summer these gains are smaller, for higher load ranges 
when the temperature of the cooling water exceeds 30°C it 
is advisable even to increase the flow of the cooling water 
with respect to the nominal mass flow rate, which can be 
done using the inlet guide vanes. 

Similar gains can be expected in case of steam turbines 
working in marine power plants. When the ship sails in the 
moderate zone and at reduced driving power we can expect 
the reduction in power plant fuel consumption amounting 
approximately to 240 kg per 24 hours. At the same time, the 
control making use of adjustable guide vane blades provides 
good opportunities for increasing the cooling water mass 
flow rate to obtain deeper vacuum in the condenser in the 
tropical zone. 


NOMENCLATURE 
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£e 
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efficiency of turbine Low Pressure (LP) part [-], 
— steam mass flow rate [kg/s], 

— cooling water mass flow rate [kg/s], 

— pump outlet pressure [m of H,O or MPa], 

— power of power unit [kW], 

change of power of power unit [kW], 

— volumetric flow rate [m/s], 

— condenser pressure [kPa], 

— power of pump motor [kW], 

— inlet temperature of cooling water [?C], 

— specific volume of condensing steam [m?/kg]. 
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Influence of hydraulic oil viscosity 
on the volumetric losses in a variable 
capacity piston pump 


Jan Koralewski, M. Sc. 
Gdańsk University of Technology 


ABSTRACT 


The variable capacity piston pumps are elements of the great power and highest energy 

efficiency hydrostatic drives. They are used in the drive systems of ship equipment such 

as deck cranes, steering gears, main propulsion of smaller vessels. The laboratory and 

simulation investigations of the influence of liquid viscosity on the variable capacity 
displacement pump energy losses have not been so far performed. 

The paper presents results of the investigations of impact the hydraulic oil viscosity has 

on the volumetric losses in a piston pump operating in the full range of its capacity and 


oil pressure. 


Key words: hydrostatic drive; displacement pumps; estimation of energy losses; impact of the oil viscosity on losses 


INTRODUCTION 


A hydraulic system with variable capacity pump, as 
a structure allowing to change the motor speed, is a hydrostatic 
drive solution with the highest energy efficiency. It is used in 
the great power systems, in the situations of prolonged system 
operation, wherever the energy saving is preferred even with 
more expensive investment and greater operation requirements. 
Examples of the ship applications are the drive and control 
solutions of a deck crane, steering gear and main propulsion 
of smaller vessels. 

It is important to know the transmission energy efficiency 
in the nominal conditions but also in the whole range of the 
operating conditions (hydraulic motor load and speed, hydraulic 
oil viscosity), particularly in the most often occurring or most 
prolonged operating conditions. 

A hydraulic drive system designer or user has at his disposal, 
provided only by some manufacturers, results of the energy 
efficiency tests of machines in the systems, tests performed 
with a selected oil viscosity. The efficiency of a hydraulic motor 
and its driving pump as elements of a hydrostatic transmission 
system should be determined as a function of the motor shaft 
speed and torque. 

There is no tool allowing to perform full energy analysis 
of the hydrostatic transmission as a whole composed of any 
different set of machines. The transmission system efficiency 
should be presented as dependent on the hydraulic motor 
speed and load, with a possibility of evaluating the impact of 
volumetric, pressure and mechanical losses, different in various 
types and sizes ofthe machines used an also evaluation of the 
impact of pressure losses in the system conduits. 


All those losses are also a function of current motor 
operating parameters and of the oil viscosity changing during 
the system operation. 

So widely understood simulation investigations require 
a suitable model of the variable capacity pump losses and 
energy efficiency, and also a model of the system efficiency with 
such pump. In order to verify the models, it is necessary to curry 
out carefully prepared laboratory investigations of the pump, 
hydraulic motor and the whole system. Such investigations, 
with the constant recommended oil viscosity v, = 35mm's" 
were performed by M. Czynski [1]. 

The laboratory and simulation investigations of 
the impact of liquid viscosity on the variable capacity 
displacement pump energy losses have not been 
performed yet. 


MATHEMATICAL MODELS OF THE 
DISPLACEMENT MACHINES LOSSES AND 
ENERGY EFFICIENCY 


Assessment of the capability of energy savings in the 
hydrostatic drive system operation requires the system losses 
to be defined. 

The simulation determination of the system energy 
efficiency may be used for the purpose in the system design and 
operation process [2]. The following factors in the simulation 
model should be taken into account: 

- the hydraulic motor speed control system structure, 
- energy losses in the system elements, 

- the pump driving motor speed decrease, 

- the system control element characteristics, 
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- load and speed of the controlled hydraulic motor, 
- hydraulic oil viscosity, changing in the system operation 
process due to change of the oil temperature. 


In order to make the transmission system efficiency 

determination method easily applicable, it is necessary to: 

1) use the computer programs for the mathematical models, 
allowing to analyse the hydraulic system efficiency as 
a function of the decisive parameters (hydraulic motor 
speed coefficient o, and load coefficient My, ratio v/v, of 
the hydraulic oil viscosity v to the reference viscosity v,), 

2) determine the values of energy loss coefficients for the 
pump, rotational hydraulic motor or hydraulic cylinder. 
Those coefficients should be clearly defined and precisely 
determined for a given displacement machine. 


The mathematical model of the displacement machine 
losses allowing to fulfil the conditions given in points 1 and 2 
above should take into account: 

a) the form and simplicity of the description, deciding of the 
possible use of that description in the system efficiency 
model, with maintaining the system efficiency precise 
assessment, 

b) description ofthe displacement machines volumetric losses, 
allowing to evaluate the impact of hydraulic oil kinematic 
viscosity changing with oil temperature, 

C) separate treatment of the mechanical and pressure losses 
in the machine. These losses increase the required torque 
on the machine (pump) shaft but the losses are of different 
character and depend on the same parameter (viscosity 
coefficient v/v.) but in a different way. 


It is necessary to perform the laboratory and simulation 
investigations in the displacement machine real operating 
conditions. The investigations should allow to verify the 
proposed models of: 

- machine volumetric losses, 

- machine pressure losses, 

- machine mechanical losses, 

in the full range of working pressure up to nominal pressure p, 
in the wide range of pump capacity up to theoretical capacity 
Q,, and in wide range of the hydraulic oil kinematic viscosity 
v, and also to determine the k, coefficients of specific losses. 


MODEL OF THE DISPLACEMENT PUMP 
VOLUMETRIC LOSSES 


Volumetric losses require an increase of the pump 
geometrical capacity, are connected first of all with the working 
liquid leaks through slots between displacement elements and 
the working chamber walls, distributor (if it exists) elements 
and are also effect of the liquid compressibility, change of the 
pump working volume and change of the slot height due to 
changes of pressure and temperature. 

The model of volumetric losses presented by Z. Paszota in 
[2, 3, 4] meets the requirements given in chapter 2. The author 
assumes the conditions and simplifications of the impact of 
certain factors on those losses and that impact is reflected in a 
coefficient and in power exponents describing the dependence 
of losses on Ap,, and v. 

The theoretical pump working volume q, (theoretical 
capacity q,, for one pump shaft revolution) — the geometrical 
difference between the maximum and minimum volume of 
working chambers — is a characteristic value of a pump. It is 
determined at the pressure value p,,, = 0 in the pump working 
chambers during their filling and at the increase Ap, = 0 of the 
indicated pressure in the chambers. 
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Flow path 


Suction conduit | Discharge conduit 


Working chambers 


Fig. 1. Approximate changes of the working fluid pressure 
along the flow in the pump 


Under the pressure and temperature, the geometric working 
volume pe of the pump changes slightly compared with q,,. It is 
assumed (in order to simplify the description of volumetric loss 
intensity Q, in the pump) that the theoretical pump working 
volume q,,is constant and equal to the geometrical working 
volume q,, determined at the working liquid temperature 
corresponding to the recommended kinematic viscosity 
v, 735 mm's*: 


Ppii70 
Apt = Apg | APri=0 (1) 


Vn 


and the change of geometric working volume d during the 
system operation will be taken into account in the values of 
loss coefficients in the pump. 

The theoretical capacity Q, of a constant capacity pump is 
given by the formula: 


Qh = Ippo (2) 
where: 


n = 


" the rotational shaft speed of an unloaded pump 


(Ap,, = 0). 


The intensity Q,, of the pump volumetric losses is 
described by a simulation model: 
-0.8 


q V 
Qp, = kpss = Appi a 


n n n 


k — dimensionless constant of volumetric losses in the 
pump, determined experimentally at the reference 
viscosity v, = 35 mm's', 

theoretical working volume of a constant capacity 
pump, 

reference mass density of the working medium 
(hydraulic oil) determined at the temperature 
corresponding to the kinematic viscosity v, and 
pressure p = 0 (atmospheric pressure), 

indicated pressure increase in the pump working 
chambers, 

kinematic viscosity of the working medium 
(hydraulic oil) used for calculation ofthe volumetric 
losses Q, and determined at the pump inlet, 
reference kinematic viscosity of the working 
medium (hydraulic oil) v, = 35mm’s", determined 
at the pressure p = 0 (atmospheric pressure), 


(v/v °- approximate description of the impact of 
liquid viscosity v on the volumetric losses in 
a displacement rotational machine. 


The value -0.8 of the exponent takes into account first of 
all two types of volumetric losses in the pump — dominating 
leaks of the laminar flow, proportional to (v/v,)' and leaks of 
the not fully developed turbulent flow, proportional to (v/v. ) ^^. 
Therefore, the -0.8 exponent may be replaced by a different 
value in a more precise description of the intensity Q,, of 
volumetric losses in specific pump. 

The value of that exponent must be determined 
experimentally for each type of a displacement pump. 


The pump capacity Q,is described by the expression: 


n n n 


-0.8 
= dp M 

Qp = qrp — Kpss —7— AP pj "i (4) 

where: the speed n, lower or equal to n,,, depends on the 

characteristic of the pump driving motor (n, decreases when 


the torque M, required by the pump increases). 


Coefficient k, of the volumetric losses Q,,, determined 
during one shaft revolution of a constant or variable capacity 
pump, at the pressure increase Ap,, equal to the hydraulic 
system nominal pressure Ap,, = p, and at the viscosity v,, the 
losses related to the pump theoretical working volume q,, is 
described by the formula: 


Q Opt 
App;— 
Pv Ppi-Pn 
oo (5) 
k, =—a—_ — 

dpt dpt 

Np APpi= Pn 
Vn 


The relation between coefficient k, and the constant value 
kps of the volumetric losses in the pump is the following: 


k, = kpas —— HÀ — — (6) 


Apt 
PrV Np Appi=Pn 


Vn 


The relation between coefficient kiy calculated at the oil 
viscosity v changing during the drive system operation and 
coefficient k, is the following: 


—0.8 
k= : (7) 
M 


n 


Itresults from the presence of two types of volumetric losses 
in the pump: dominating leaks ofthe laminar flow and leaks of 
not fully developed turbulent flow. 

Using the volumetric loss coefficient k, the following 
formula for the intensity Q,, of volumetric losses is obtained: 
1 -0.8 
pt 
Appi =Pn 


Vn 


and the pump capacity Q, formula: 


R 1 —0.8 

dp, 

tn e: M m 
p 


v ü V 


n 


App 


n 


= (8) 


Q,, = kqyn, 
V 


n 


Q,-qyny-kgqgan, 


The use of coefficient k, of volumetric losses for description 
of the relation of Q,, intensity to the indicated increase Ap,, 


of pressure in the pump working chambers allows to describe 
that relation by an exponential function with the exponent not 
necessarily equal to 1. The use ofk,,.. constant to the description 
of Q,, required an assumed proportionality of Q,, to App; 

The expressions describing the variable capacity pump 
capacity takes the form: 


-08 
q V 
Qp= bqunp- ks; —R- Ap», =| (10) 


n n n 


A 1 —0.8 

pt SLM 

imet | 7 =| 
Vn Pa Va 


The expressions (10, 11) assume, that change of the pump 
capacity setting b, (change of the pump capacity) does not 
influence the intensity of pump volumetric losses Q,,. 

In the expressions (9, 11) the value 1 of exponent describing 
the impact of the Ap,/p, ratio and also the value -0.8 of exponent 
describing the impact of the v/v, ratio on the intensity Q,, 
of pump volumetric losses should take into account all the 
factors influencing the volumetric losses (character of the 
flow in slots, change of the slots cross-section with pressure 
and temperature, liquid compressibility, change of the liquid 
viscosity in slots etc.). 

The value 1 of exponent describing the 
impact of the Ap,/p, ratio and also the value -0.8 
of exponent describing the impact of the v/v, ratio on the 
intensity Q,. of pump volumetric losses must be verified 
experimentally for each pump type. 


RESULTS OF THE LABORATORY 
INVESTIGATIONS 


Laboratory investigations of an axial piston variable 
displacement pump of bent axis design (BOSCH REXROTH 
A'IVS8RD type) were carried out on a test stand in the Chair 
of Hydraulics and Pneumatics of the Gdansk University of 
Technology Mechanical Engineering Faculty. 


or: 


Q,-b,qyn,- k,qp.Np 


Fig. 2. Axial piston variable displacement pump of bent axis design 
(BOSCH REXROTH A7V58RD type) 
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The investigations were performed with: - 7 values of pump capacity coefficient: 


- 8 temperatures of hydraulic oil (oil kinematic viscosity v): b, = 0.227; b, = 0.361; b, = 0.493; 
v —20*C (v= 120.40 mm’s"), v = 24°C (v= 91.16 mm’s"), b, = 0.623; b, = 0.752; b, = 0.880; b, = 1. 
v = 30°C (v= 65.37 mm’s"), v = 36°C (v = 47.05 mm’s’"), 
v = 43°C (v = 34.68 mm’s"), v = 50°C (v = 26.41 mm’s’"), The selected method of determination of the geometrical 
v = 60°C (v = 18.77 mm’s"), v = 68°C (v = 14.53 mns), (variable) capacity q,,, per one pump shaft revolution and 
- 8 values of the increase Ap, of pump pressure: theoretical capacity q,, per one shaft revolution was based 
Ap, = 1.6 MPa, Ap, = 3.2 MPa, Ap, = 6.3 MPa, on the extrapolation of linear functions q, = f(Ap,,), in the 
Ap, = 10 MPa, Ap, = 16 MPa, Ap, = 20 MPa, range of small increases Ap,, of pressure in the pump working 
Ap, = 25 MPa, Ap, = 32 MPa chambers (Fig. 3). 
bp=0.227 Oil viscosity ratio v/v, v.-35mm's! — bp-0.493 Oil viscosity ratio v/v, va=35mm s" 
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Fig. 3. Determination of the geometrical (variable) capacity q PN (ey = b,q,) per one shaft revolution and of the value of the pump capacity coefficient b, 
from the relation of the pump capacity q, per one shaft revolution to the indicated increase Ap, of pressure in the pump working chambers at different values 
of the oil viscosity ratio v/v; examples for four choices of different pump capacity settings b, = 0.227 + 1 
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Description of q, = f(Ap,,), with those linear functions in 
the range of small pressure increases Ap,,, allowed to determine 
poy (Ap) with the accuracy of an order of 1 per mille (0.001). 
Approximation, instead with a linear function in the whole 


range of the increase Ap,, of pressure (up to 32 MPa) or a 
second degree polynomial or an exponential function in the 
whole or a small range of Ap,, allowed to determine ds. (CH) 
with much less accuracy. 
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Fig. 4. Relation of the intensity Q,, of volumetric losses to the indicated increase Ap of pressure in the pump with different values of the oil viscosity ratio 
v/v ; examples for four chosen values of the pump capacity coefficient b, = 0.227 + 1. In the range up to Ap,, = 16 MPa, the Q,, intensity is most precisely 
described by the linear functions, in the Ap,, = 16 + 32 MPa range, the Q,, intensity is described by exponential functions 
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The figure 4 demonstrates a complex impact on Q,, of the 
flow character in the slots and of the changes of slot cross- 
section and the hydraulic oil viscosity under the influence of 
pressure and temperature. 
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with the measurement results at Ap,, = p,. A consequence of 
the choice of such functions is worse conformity with the 
measurement results at lower Ap,, values. 
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Fig. 5. Relation of the intensity Q „, of volumetric losses to the indicated increase Ap,,of pressure in the pump, described by exponential functions in the 
whole range of pressure; examples for four chosen pump capacity coefficients b, = 0.227 + 1 
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Fig. 8. Relation of the intensity Q,, of volumetric losses to the pump capacity coefficient b, at different values of indicated increase Ap, of pressure in the 
pump working chambers and at different v/v, = 0.42 + 3.47 ratios of oil viscosity 
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Fig. 9. Relation of the intensity Q,, of volumetric losses to the oil viscosity ratio v/v, at different values of indicated increase Ap, 
of pressure in the pump working chambers and at different values of the pump capacity coefficient b, = 0.227 + 1 
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Fig. 11. Value of the a, exponent (in the exponential function describing 
the relation of the intensity Q,, of volumetric losses to the oil viscosity ratio 
v/v.) at changing pump capacity coefficient b, and for different values of the 
indicated increase Ap, of pressure in the pump 


Values of the a exponent, in the exponential function 
ERA apy H $ B 
describing the relation ofthe intensity Q,, of volumetric losses 
to the indicated increase Ap,, of pressure in the pump, presented 
in figures 6 and 7, are within the 0.91 « a, < 1.31 range. This 
range is limited to 0.91 < a, < 0.96 at the oil viscosity ratio 
v/v, — 3.47. 

"The exponent values a , < 1 obtained at the highest oil 
viscosity allow to conclude that in the whole range of viscosity 
v the flow in slots is of a not fully developed turbulent character 
(with increasing turbulence at decreasing viscosity). Increase of 
the a, exponent above | at decreasing viscosity v indicates an 
impact of the slot increase upon the intensity Q, of volumetric 
losses as an effect of temperature increase. 

The pump capacity coefficient b, has practically no impact 
on the intensity Q,, of volumetric losses in the pump working 
chambers (Fig. 8). 

Values of the a , exponent in the exponential function 
describing the relation of the intensity Q,, of volumetric losses 
to the oil viscosity ratio v/v „ presented in figures 9, 10 and 11, 
in the a = -0.20 + -0.35 range indicate the domination of not 
fully developed turbulent flow over laminar flow in the pump 
slots. 


VERIFICATION OF THE MATHEMATICAL 
MODEL OF PUMP VOLUMETRIC LOSSES 


In order to verify the mathematical model described by 
formula (8), it was replaced by a mathematical expression taking 
into account the relations, obtained during the investigations, 
of the intensity Q,, of pump volumetric losses to the indicated 
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Fig. 12. Determination of the a, exponent in the mathematical model 
describing the relation of the intensity Q,, of volumetric losses to the 
indicated increase Ap, of pressure in the pump working chambers; pump 
capacity coefficient b, = 1, the oil viscosity ratio v/v, = 1. From formula (5), 
the value of coefficient k, — 0.065 of volumetric losses is determined. 
The value of a= 0.97 exponent is obtained 


increase Ap,, of pressure in the pump working chambers (to the 
Ap,,/P, ratio) and also to the oil viscosity ratio v/v: 
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Q= kqyn, v, SENI (12) 
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Assumed was (Fig. 12) the value of exponent a, = 0.97 in 
the formula (12) determined with the pump capacity coefficient 
b, — ], the oil viscosity ratio v/v, — ] and the coefficient of 
volumetric losses k, = 0.065 calculated from formula (5). 

Assumed was (Fig. 13) the value of exponent a, = -0,30 in 
the formula (12) determined with the pump capacity coefficient 
b, = 1, Ap,/p, = 1 ratio and the calculated coefficient of 
volumetric losses k, = 0.065. 

The obtained values ofthe coefficient k, — 0.065 of intensity 
Q,, of volumetric losses, exponent a, = 0.97 of the relation 
of intensity Q, of the volumetric losses to the ratio Ap,,/p, of 
the pressure increase, exponent a, = -0.30 of the relation of 
intensity Q, of the volumetric losses to the v/v, ratio of oil 
viscosity, have made it possible to present the mathematical 
model of the intensity Q, of pump volumetric losses in the 
form: 
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Model (13) describes precisely the intensity Q, of 
volumetric losses in the nominal conditions of the pump 
operation, i.e. at the pump capacity coefficient b, — 1, the 
pressure increase ratio Ap,/p, = 1 and the oil viscosity ratio 
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Fig. 13. Determination of the a. exponent in the mathematical model 

describing the relation of the intensity Q „, of volumetric losses to the oil 
viscosity ratio v/v ; pump capacity coefficient b, = 1, indicated increase 

Ap, =P, = 32 ? MPa of pressure in the pump 9 working chambers, 

coefficient k, = 0.065 of volumetric losses. The value a, = -0.30 

of the exponent is obtained 
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Fig. 14. Comparison of the intensity Q „, of volumetric losses described by the mathematical model (13) with the laboratory investigation results and the 
absolute difference between the mathematical model values and the laboratory investigation values; assumed were: the coefficient k, = 0.065 of volumetric 
losses, exponent 2,7 0.97, exponent a, = -0.30; examples for four chosen values of the pump capacity coefficient b, = 0.227 + 1 
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v/v, = 1. At the same time this model is a simulation formula 
describing the change of intensity Q,, of volumetric losses 
with the change of the pressure increase ratio Ap, /p, and the 
oil viscosity ratio v/v, (the change of pump capacity coefficient 
b, has practically no impact on the change of intensity Q,, of 
volumetric losses). 

Figure 14 presents a comparison of the intensity Q,, of 
volumetric losses described by the mathematical model (13) 
with the results of laboratory investigations, supplemented 
by the information about the absolute difference between 
the values from the mathematical model (13) and results of 
laboratory investigations. Examples are given for four selected 
values of the pump capacity coefficient b, 

Differences between the simulation and experimental values 
of the intensity Q,, of volumetric losses, determined in the 
whole range of the pressure increase ratio Ap, /p,, oil viscosity 
ratio v/v. and pump capacity coefficient b, are mainly caused 
by the change of a, exponent describing the relation of the 
intensity Q,, of volumetric losses to the pressure increase ratio 
Ap,/p, in the situation of using in the mathematical model the 
value a, = 0.97 determined at b, = 1 and v/v, = 1. 


CONCLUSIONS 


1. The purpose of the investigations was experimental 
verification of the mathematical model (8) [2, 3, 4], 
describing the volumetric losses in a variable capacity 
displacement pump used in hydrostatic transmissions. Model 
(8) allows to describe the losses and the energy efficiency 
of the pump and hydrostatic drive as a function ofthe drive 
speed and load and also the hydraulic oil viscosity. 

2. Model (8) allows a simple and precise determination ofthe 
pump volumetric losses by determining the coefficient k, 
of volumetric losses (5) in the nominal conditions of pump 
operation — at Ap, = p, bp = 1, v/v, = 1. 

3. Model (8) allows also to determine the impact of the ratio 
Ap,,/P, of pressure increase, the hydraulic oil viscosity ratio 
v/v, on the intensity Q,, of volumetric losses in the whole 
range of the pump capacity coefficient b, 

4. The investigations were carried out with an axial piston 
variable displacement pump of bent axis design, commonly 
used in hydrostatic transmissions. 

5. In order to verify the mathematical model (8), it was 
replaced by formula (12) for investigating the exponent a, 
in the expression Qy” (Ap, /p,)"»* and exponent a, in the 
expression Q,, ~ (viv, J^. 

6. The chosen method of determining the pump geometrical 
working volume s. and theoretical working volume q,, 
was based on extrapolation of linear functions q, = f(Ap,,) 
within the range of small pressure increases Ap, in the 
working chambers. This allows to determine Wi (dp) with 
the accuracy of an order of | per mille (0.001). 

7. Acomplex impact on the intensity Q,, of volumetric losses 
was found of the character of flow in slots as well as the 
impact on Q,, of the change of slot cross-sections and 
hydraulic oil viscosity v due to the change of pressure and 
temperature. Up to Ap,, — 16 MPa, the intensity Q,, was best 
described with linear functions, in the Ap,, = 16 + 32 MPa 
range the intensity Q,, is described by exponential functions. 
For description ofthe Q, to Ap,, relation in the whole range 
of the pressure increase, the exponential functions giving 
the best agreement with the mesurement results in the Ap,, 
=p, area were chosen. 

8. The values of exponent a, in the expression Q,, ~ (App, )' 
are in the 0.91 « a, < 1.31 range narrowing to the 
0.91 <a, < 0.96 range ' at the oil viscosity ratio v/v. — 3.47. 


9. It has been found out that the pump capacity coefficient 
b, has practically no impact on the intensity Q,, of pump 
volumetric losses. 

10. The values of exponent a, in the expression Q,, ~ (v/v ^w 
are in the -0.35 «a. < -0. 20 range and show the domination 
of not fully developed turbulent flow over the laminar 
flow in the pump slots in the whole range of investigation 
parameters. 

11. The value k, = 0.065 of the coefficient of the volumetric 
losses in the pump working chambers was calculated at 
the pump capacity coefficient b, = 1, the pressure increase 
ratio Ap,/p, = 1 and the oil viscosity ratio v/v, = 1. The so 
determined value of k, coefficient allows the quantitative 
and qualitative evaluation of the pump volumetric losses. 

12. The value a = = 0.97 of the exponent in the expression 
Qu ~ (AP, Jp, )* was determined at the pump capacity 
coefficient b = = | and the oil viscosity ratio v/v, = 1. 

13. The value a. — -0.30 of the exponent in the 'expression 
Q,- (v/v yw was determined at the pump capacity coefficient 
b, = 1 and the pressure increase ratio Ap, /p, = 1. 

14. In effect, the mathematical model of volumetric losses in 


the investigated pump takes the form (13): 
0.97 -0.30 
V 


Opt : 
Q,-0065q,n,. (AP 
Appi =Pn Pa Va 
15. Intensity Q,, of pump volumetric losses described by the 
mathematical model (13) was compared with the results 
of laboratory investigations. The absolute difference 
between the values from the model and from the laboratory 
experiment did not exceed: at b, = 1 — +2 + -6 cms", at 
b, = 0.227 — -8 + -llem'?s' compared with the nominal 
operation value of Q,, = 94 cms". 
16.It must be underlined, that in the assumed conditions 
of determination of the k, coefficient (conclusion 11), 
the difference between the value of the intensity Q, 
of volumetric losses from the model and the results of 
laboratory investigations in nominal conditions (b, — 1 
Ap,/p, = 1, v/v, = 1) is equal to zero. 
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ABSTRACT 


The ability to determine the bearing life time is one of the main purposes in maintenance of rotating 
machineries. Because of reliability, cost and productivity, the bearing life time prognostic is important. In 
this paper, a stiffness-based prognostic model for bearing systems is discussed. According to presumed model 
of bearing and fundamental of damage mechanics, damage curve approach is used to relate stiffness of 
vibratory system and bearing running life. Furthermore, using the relation between acceleration amplitude 
at natural frequency and stiffness, final relation between acceleration amplitude at natural frequency and 
running life time according to damage curve approach can be established and the final running time is 
predicted. Experiments have been performed on self alignment bearing under failures on inner race and 
outer race to calibrate and to validate the proposed model. The comparison between model-calculated 
data and experimental results indicates that this model can be used effectively to predict the failure lifetime 
and the remaining life of a bearing system. 


Keywords: Rolling bearing, Prognostic, Damage mechanics, Damage curve approach 


INTRODUCTION 


Nowadays, rotating machineries are widely used in 
mechanical equipments and different industries. The failure 
of a rolling element of bearing is one of the primary causes 
of breakdown in rotating machineries. Consequently, 
development of reliable techniques for the prediction of rotating 
element failure is more important to facilitate preventive 
maintenance. 

The ability to determine the bearing life time is one of 
the main purposes in maintenance and condition monitoring 
of rotating machineries. If remaining life time of bearing be 
predicted then reliability of bearing will be improved and 
cost of maintenance will be reduced, so bearing life time 
prognostics is one of key points of condition monitoring of 
rotating machineries. 

Maintenance of bearing can be categorized in two main 
group of diagnostic and prognostic. Currently, there are three 
methods for bearing diagnostic. These methods are vibration 
analysis [1, 2, 3], statistical study [4, 5] and neural network 
[6, 7]. The researches in the field of diagnostic of failures of 
bearings are focused on determination of failures of bearing 
and do not predict the running life of bearing. 

Diagnostic do not estimate the bearing life while prognostic 
is focused on prediction of bearing life time. The studies of 
prognostic of remaining life of bearing is categorized in two 
approaches: Direct approach of fracture mechanics [8, 9, 10], 
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and indirect approach which is a statistical analysis of vibration 
signals [11]. 

The primary researches to estimate the remaining life of 
bearing was carried out in 90 decay. Li et al., 1999, used the 
vibration analysis and fracture mechanic to produce a recursive 
algorithm to predict the remaining life of bearing. These 
research shows that fracture mechanics can be useful to predict 
the bearing life. Qiu et al., 2002, used damage mechanics to 
predict the bearing life. In this research, three approaches of 
damage mechanics [Linear Damage Rule (LDR), Damage 
Curve Approach (DCA) and Double Linear Damage Rule 
(DLDR)] for trend the variation of acceleration amplitude 
on natural frequency of system were compared. However, 
the deterministic formula for estimation of remaining life of 
bearing did not explain. Due to previous studies in the fields 
of diagnostic and prognostic of bearing remaining life, this is 
obvious that most of researches are focused on diagnostic of 
defects of bearing and prognostic of bearing life time is a new 
branch in the field of condition monitoring and maintenance 
of bearings. 

As a result, there is still a lack of reliable prognostic 
methods to predict the remaining life of berating bearing. 
Reliable and physically meaningful estimation of bearing life 
under various operating conditions remains a challenging and 
important aspect for maintenance optimization and catastrophe 
avoidance. In this paper, a stiffness-based prognostics model 
for a bearing system based on vibration response analysis and 


damage mechanics is developed. The effectiveness of the model 
has been studied by comparison between estimated life due to 
model and real running life of experiment, and accuracy of the 
model is also investigated. 


BEARING DYNAMIC MODEL 


A bearing system can be introduced with a self-excitation 
dynamic system. In this system, two types of loads exist. One 
type is external loads which affect bearing because of working 
situation and the other loads are due to defect of bearing such 
as surface defects, unbalancing, etc. According to the studies of 
Qiu et al., 2002, bearing system can be simplified as single-dof 
systems as shown is figure 1. 

In this model, the stiffness and damping of system have 
two parts: undamaged part and damage part. Undamaged part 
is the part of bearing which is remain safe during test running 
life and damage part is the part of bearing which is finally cause 
bearing break down. According to the fracture mechanics, high 
cycle fatigue is the main reason of bearing failure and during 
bearing life the two phase of crack (initiation and propagation) 
are happened. Damage part of bearing in this model is the 
presumed crack which finally leads to failure and its stiffness 
reduced during the working life. The total stiffness and total 
damping of bearing system 1s summation of stiffness and 
damping of damaged part and undamaged part, respectively, 

K *K, & C-C, * C, In this equation u denotes undamaged 
part and d denotes damaged part. 


Fig. I. The dynamic model of bearing system 


Because of rotating nature of bearings the forces which are 
affected on bearing systems are periodic, and we call them F(t). 
Expanding F(t) by Fourier transform, the acceleration response 
of system can be obtained as: 


—Q?F;sin (Q;t + aj) 


x= —————— (1) 
i |(k — MQF)? + (C0)? 
where: 
Q, — the frequency of the ith harmonic excitation 
F,á — the phase angle between F, and x. 


The amplitude of acceleration response at the natural 
frequency is: 
OnFn 


a - (2) 


Presume proportional damping for system, C = aM + bK 
and considering the nature of presumed system, such as no 
mass enters and no mass comes out of system and the mass 
of system is constant, then a = 0. Consequently, damping of 
system is only a proportion of the system stiffness. With this 
presumption, equation (2) can be written as: 


C 
A, == (3) 
” K 
where: 
c=F/ (b/M) — a parameter which is depend on characteristics 
of system. 
DAMAGE MECHANICS MODEL OF 
BEARING 


Under the definition ofthe reference volume element (RVE) 
in damage mechanics of material[12], and previous studies in 
this field (Qiu et al., 2002), a damage factor can be generally 
defined as: 


Ka 
D-d (1 - =") (4) 
Kao 

Where d is scaling factor [d = K,,/(K,, — Kag and K, is the 
stiffness of damage part which is reduce during the working 
life of bearing and initial and final value of it are K,, and K p 
respectively. Based on damage mechanics, Damage curve 
approach can be used to relate the running life of bearing with 
damage factor as below 


N q 
D= (x) (5) 


where: 

N - therunning cycles, 

N,- the failure lifetime in cycles, 

q — amaterial and structure-dependent factor. 


From equations (4) and (5), it can be shown that: 


Ka 1/N\9 
ze 717a) © 
Kao d AN; 
On the other hand, the equation (3) can be written as 
below: 
| |. Co ( Ka + Ku ) T 
Aa i Ca Kao + Ky 


Where A, and A, are the amplitudes of acceleration 
response at the natural frequency under damaged and initial 
conditions, respectively, and c, and c, stand for coefficient c in 
equation (3) under damaged and initial conditions, respectively. 
Combination of equations (6) and (7), the final equation for 
estimation of the failure life is obtained as below: 


(S = eha - TAn x5) (x) ) = 
Nd (8) 
=a (: — By) 


Where á, à and q are the coefficients depending on the 
operating condition, materials and structure of the system. 
Subscript d denotes the damaged condition, subscript 0 denotes 
the initial condition and subscript f denotes the final life time 
of bearing. For bearings, the failure time represents the time 
when the vibration signals tend to be infinite. In equation (8), 


â and q are constant for each bearing then with substitution 
Nfe = a equation (8) can be written as below: 
vB 
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It is obvious that if vibration signal tends to be infinite 
then N will tend to N,,. According to equation (9), obviously 
if unknown parameters, á and q be determined, and amplitude 
of acceleration response at any time of useful life of bearing 
be measured, then the lifetime of bearing system can be 
calculated during running life of bearing. To determine the 
unknown parameters, Recursive Least Square (RLS) parameter 
estimation algorithm can be used. Following, an experimental 
test rig is designed to validate the theoretical equation. Using 
RLS algorithm and experimental data, first the unknown 
parameters are identified and then accuracy of estimated life 
of bearing is investigates. 


EXPERIMENTAL TEST RIG 


For prognosis bearing life time an experiment test rig is 
used. The experimental test rig has two parts, electrical and 
mechanical parts. The electrical part is data acquisition system 
which is used to acquire the bearing vibration signals. The 
mechanical part of experimental test rig 1s shown in Fig. 2. 
Two supporting bearings, shaft, electromotor and main bearing 
test are segments of mechanical part of test rig. The shaft was 
driven by an electromotor with a speed controller. 

Supporting bearings are self-alignment spherical roller 
bearing 22311EXK. These bearings can tolerate the fatigue 
stresses during fatigue test of bearing and would not fail 
during the experiment. The main bearing which is studied is a 
self-alignment ball bearing 1206C3. Defects on inner race and 
defects on outer race are main defects of this type of bearing, 
so both of these defects are studied. The operating conditions 
of the experiment are: shaft speed 1500 rpm and radial load 
is equal to 9 kN. 


! este d bearing 


EE pos — c: 


— i; 


Two supporting Load Indicator 


Fig. 2. Experiment test rig 


The radial load which is applied to the tested bearing is 
measured by a load-cell (the load-cell in connected to load 
indicator which shown the value of load). The manufacturer of 
Load-cell and load indicator is Bongshin Company, load-cell 
model is DBBP and load indicator model is BS-7220. In Fig. 3 
the loading configuration is shown. As shown in this figure, the 
loading configuration consists of a lower metal plate (thickness 
= 20 mm), upper metal plate (thickness = 20 mm), load-cell, 
4 machine bolts (M24) and tested bearing housing. In Fig. 3 
this parts are shown. 

According to Fig. 4, by screwing the nut of 4 machine bolts 
the lower plate goes down and consequently the load-cell which 
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is placed between two plates will be tensioned. The tension 
force leads to a radial force on bearing by bearing housing. In 
Fig. 4 the radial force is shown. The vibration signal will be 
measured by piezoelectric sensor. The manufacturer of sensor 
is ENDEVCO Company and the model of sensor is 2235C. 
The signal amplifier is made by Brüel & Kjær and its model 
is 2525. Moreover, AD convertor is made by ADVANTECH 
and its model is USB4711. In this experiment acceleration of 
tested Deanne | is eee 


lachine Holts 


Lowergplate 


Fig. 4. Loading on tested bearing 
EXPERIMENTAL RESULT ANALYSIS 


As mentioned earlier, the equation for estimation of 
remaining life of bearing has two unknown parameters which 
must be defined with RLS algorithm. For this reason four 
primary tests are carried out to calculate these parameters. In 
test 1 and 2, defect is on outer race and in test 3 and 4 the defect 
is on inner race of the bearing. After primary tests, the unknown 
parameters are determined. Equation (9) is rewritten for both 
defects on the outer race and t on the inner race, respectively, 


as follows: 
A ( N ) 
-125 1-(-) 
Oe Nfe 


(10) 


Anoy N \° 
= e mus (11) 
(72) =09(1- (52) ) 


Now, two tests are carried out to evaluate the precision of 
equations (10) and (11). In test 5, defect is on the outer race 
and in test 6, defect is on the inner race. Equations (10) and 
(11) can be written as below: 


(12) 


(13) 


In tests 5 and 6, if the acceleration amplitude at natural 
frequency of each cycle of running life of bearing be calculated 
and consequently ratio of calculated amplitude acceleration to 
the initial amplitude ((A_,/A, )?^) determined, then the remaining 
life of bearing can be estimated. 

Before estimating the life of bearing in tests 5 and 6, the 
rating life of bearing evaluated. The theoritical life of bearing 
calculated according to ISO 281:1990 by following equation: 


C 10/3 
Lio = (5) (14) 


The C parameter for tested bearing 1s equal 15600 N and as 
mentioned previously the P value is equal 9000 N; consequently, 
theL, , equals to 6.25* 10° revolutions. The obtained estimations 
of life of bearings are shown in figures (5) and (6). In figures 
(7) and (8), error percentage of the estimated life in tests 5 and 
6 are shown. 

As shown in figures (7) and (8), the final life of bearing 
estimated precisely during running of bearing. In test 5, 
during 31 percent of last running life of bearing the error of 
estimated life of bearing is below 15 percent. Moreover, in test 
6, during 37 percent of last running life of bearing the error 
of the estimated life of bearing is below 15 percent. The final 
estimation is more precise, error percentage of final estimation 
of life of test 5 is 0.36 and for test 6 the final error is 1.17 
percent. Consequently, damage curve approach can be used 
for bearing prognostic in an online manner. 
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Fig. 5. Estimated life of the ball bearing in Test 5 (defect in outer race) 


6 Estimated life of Test 6 
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Fig. 6. Estimated life of the ball bearing in Test 6 (defect in inner race) 
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Fig. 7. Error percentage of estimated life in Test 5 


Error percentage of Test 6 
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Fig. 8. Error percentage of estimated life in Test 6 
CONCLUSION 


This paper is used the damage curve approach to estimate 
the remaining life of bearing. Using damage mechanics and 
single-dof model for bearing system, the appropriate equation 
for estimation offinal life of bearing is developed. Experiments 
are carried out to calculate the unknown parameters of 
theoretical equations by RLS algorithm during four primary 
tests. After primary tests, the unknown parameters determined. 
Finally, two tests are carried out to evaluate the precision of 
equations, and the final life of bearings estimated precisely 
during running of bearing (defect on inner race and defect 
on outer race). During last third of running life of bearing the 
error percent of estimated life is below 1596. Consequently, 
the damage curve approach can be used to estimate the final 
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life of bearing in online condition monitoring program. Doing 
this, firstly, using the DCA and RLS algorithm the predictive 
mathematical model must be made and then, for the other 


be 


arings we can used this mathematical model to estimate the 


remaining life of bearing precisely. 
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